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Simplified Estimators for the Normal Distribution 


When Samples Are Singly Censored or Truncated* 


A. Cuirrorp CoHEN, Jr. 
The University of Georgia 


Easily computed estimates of the mean and variance of a Normal distribution 
obtained from samples which are truncated or censored are described. The estimates 
require only a single auxiliary function which is conveniently tabulated. Worked 
examples are given for sample data subject to left and right truncation. Additional 
examples are given for left and right censoring for the two cases in which either the 
point of censoring is fixed or the number of censored items is known. 


1. INTRODUCTION 


In life testing, dosage-response studies, target analyses, biological assays, 
and in other related investigations, sample selection or observation is often 
restricted in some portions of the range of possible population values (i.e. in 
certain intervals of the random variable). Depending on the nature of these 
restrictions, resulting samples are designated as truncated or as censored. 
Samples from which some population values are entirely excluded are described 
as truncated. Those in which sample specimens whose measurements fall in the 
restricted intervals may be identified and thus counted but not otherwise ob- 
served, while the remaining sample specimens may be observed without restric- 
tion, are designated as censored. This paper is limited to maximum likelihood 
estimation of the mean and variance of the normal frequency distribution with 
probability density function 


f(a) = (0-V 2m)" exp —(x — u)’/20?, -o Sado, (1) 


when observation in samples of these types is restricted in only one ‘tail’ of the 
distribution. For the convenience of readers who might wish to explore the 
subject more thoroughly, a list of references to some of the previous work in 
restricted sampling is appended. 

Estimators derived here are simpler in form and can be more readily evaluated 
than maximum likelihood estimators previously given for the cases under con- 
sideration. They are expressed as restricted sample means and variances respec- 
tively, plus simple easily computed corrections which involve only a single 
auxiliary function of the sample terminus. Calculation of estimates accordingly 


* Sponsored by the Office of Ordnance Research, U. 8. Army. Presented to the International 
Statistical Institute in Brussels, Belgium on Sept. 4, 1958 and to the American Statistical 
Association in Chicago, Ill. on Dec. 28, 1958. 
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involves interpolation in only one table. With the exception of estimators given 
by Gupta [7], who considered singly censored samples only, previous applicable 
maximum likelihood estimators have involved two or more auxiliary functions 
and therefore interpolation in two or more separate tables. Estimators derived 
by Ipsen [10] for singly censored samples from a normal distribution also involve 
only a single auxiliary estimating function and thus interpolation in only one 
table. However, his estimators, which are based on certain moment functions 
of the restricted distribution, differ slightly from applicable maximum likelihood 
estimators. Furthermore, his tabular intervals are too wide and his entries 
contain too few significant digits for accurate interpolation. Gupta’s maximum 
likelihood estimators employ an auxiliary function which unfortunately lacks 
linearity even over short intervals of his argument. Consequently, his tabular 
intervals also are in many instances too wide for easy interpolation. Auxiliary 
functions employed here are approximately linear over moderately wide intervals 
of the arguments for both truncated and censored samples, so that accurate 
interpolation between table entries is relatively easy in both cases. Tables and 
graphs of these auxiliary functions are appended. 

Since estimators of this paper were derived by the method of maximum likeli- 
hood, for a given sample they lead to identical estimates except for possible 
errors of calculation, that might be obtained from applicable maximum likelihood 
estimators previously derived by Fisher [6], Hald [8], Halperin [9], Gupta [7], 
the author [1], and perhaps by others. The computing routine given here, however, 
is believed to be much simpler and easier to carry out. As with maximum likeli- 
hood estimators in general, those for truncated and censored samples are con- 
sistent and asymptotically efficient. They are to be recommended when sample 
sizes are at least moderately large. When estimates must be based on samples 
of size 10 or less—perhaps even on slightly larger samples, it might be preferable 
to employ linear unbiassed estimators based on order statistics as given by 
Gupta [7] in the latter part of his paper and by Sarhan and Greenberg [12]. 


2. Sminaty TRUNCATED SAMPLES 


Let x) be a known fixed value of the random variable, x, which we designate 
as a terminus or truncation point. Now consider a sample consisting of n observa- 
tions (values) of this random variable, such that for each observation (i.e. for 
each sample value), either 


(a) x > 2, in which case truncation is on the left, 


(b) x < 2, in which case truncation is on the right. 


The number of otherwise possible sample values excluded from observation as a 
consequence of this restriction is not known. 

Since probability density function (1), is symmetrical about py, truncation 
of f(x) on the right at x is equivalent to truncation of f(—2) on the left at —2 . 
It is therefore necessary to examine only one of these cases in detail, and for 
this role, truncation on the left has been selected. 
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Let F(x) designate the distribution function of x, and the probability that a 
selected value of this random variable meets the requirements for inclusion 
in a sample that is singly truncated on the left at xo is given as 1 — F(x») or in 
standard units as 1 — F(é), where 


F@) = [ o dt, with § = (%— s/o, and 


g(t) = (W2n)” exp —2/2. 


The likelihood function for a sample of the type under consideration is 


Plo ti y *** uj Hy 0) = [1 — FOV" V2n)™ exp |- ; (a; — 1/20" |. (3) 


Maximum likelihood estimating equations follow as 


(2) 


Lo — wp = Ob, 


E—p=cZ, 
s+ (€ — »)’ = o[1 + £2], 


where # and s° are the sample mean and variance respectively 


g= i and 3’ = > (x; — #)’/n, 


and where 

Zé) = o@)/[1 — F@). (5) 
The first equation of (4) follows from the second equation of (2). The last two 
results from taking logarithms of (3), differentiating with respect to u and o° 
in turn, and equating resulting derivatives to zero. The required estimators, 
A, é and the auxiliary estimator £ are to be found as simultaneous solutions of 
(4) in terms of the sample statistics, Z, x) , and s. Throughout this paper, the 
symbol (~) serves to distinguish maximum likelihood estimators from the 
parameters being estimated. 


On eliminating (¢ — y») between the last two equations of (4), we have 
=o [1-22-28] o & =8 +072 — 8). (6) 
Eliminating » between the first two equations of (4) leads to 
E—%=0(Z —£) or o = (€ — X%)/(Z — 8). (7) 
Combining (6) and (7) gives 


a =8 + (Z/(Z — d\( — 2%)’. 
Now let 


Z(é) 
a €° 


and the estimating equation for o” assumes the form 


o = 8 + WF — %)’. 


ag) = 
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To derive a corresponding equation for estimating » which does not involve 
any auxiliary function other than 0, we eliminate (Z — £) between (6) and (7) 
to obtain ¢Z = (c° — s*)/(& — xo). On combining this result with (9), we have 


oZ = W€ — 2X). (10) 


When (10) is substituted into the second equation of (4), we write the desired 
estimating equation as 


w= £— OF — x). (11) 


By eliminating o between (6) and (7), we obtain the more familiar result, 
see for example reference [1] 


[1 -— 22 -Hl/(Z — & = # /@ — x)’. (12) 


The system of estimating equations (4) may now be replaced by the equiv- 
alent system consisting of (9), (11), and (12). Let £ designate the solution of 
(12), let 6 = 0(£), and the desired estimators become 


= 8° + W% — x)’, 


(13) 
p=z— 6% —2X). 


As computational aids, tables, and a graph of @ as a function not of £, but of 
[1 — Z(Z — £)]/(Z — &)? have been provided. Since £ is the value of ¢ for which 
[l — Z(Z — &)]/(Z — &)? = s°/(& — 2x)”, we can thereby determine 6 directly 
for any given sample as that value of @ which corresponds to the sample statistic 
s’/(& — 2)°. Accordingly, the necessity for determining £ explicitly prior to 
calculating 6 is eliminated, and since @ is the only auxiliary function appearing 
in estimators (13), only the single table of that function is needed in contrast 
to the two or more tables necessary when employing estimators previously 
proposed. 

Entries of @ in Table 1 were computed from existing tables of normal curve 
areas and ordinates at equal intervals of ¢. This, of course, resulted in unequal 
intervals of the argument s”/(@ — x)”. Although equal intervals of this argument 
might be desirable, a degree of accuracy adequate for most practical applications 
can be achieved through simple linear interpolation, and the table has proven 
easy enough to use even with unequal intervals. In view of this fact, and since 
any graduation for the purpose of equalizing intervals would either result in a 
loss of significant digits or require complete recomputation of the table, it is 
being published in its present form. 

With x, , , s’, and accordingly s’/(@ — 2»)” available from the sample data, 
it is necessary only that we read 6 from the table or graph as required and 
calculate 6° and a from (13). In many applications, 6 may be read with sufficient 
accuracy from the graph of Figure 1. When more accurate values are required, 
they can be obtained by direct reading or by linear interpolation from Table 1. 
Only in rare cases should it be necessary to resort to more complicated non-linear 
interpolative procedures. 

Once 6° and @ have been computed, £ follows from (2) without the need of 
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additional tables as 


E a (xo — A)/é, 


Ve. 


where 


TaBLeE 1 
Auziliary estimation function, 0, for singly truncated samples 


8/(Z — x0) 3*/(Z — x) 


0.062 46 0.155 82 0.008 09 


-064 05 ; -156 86 -008 32 
.065 69 : .157 90 -008 56 
.067 39 3 .158 95 -008 81 
-069 16 ; -160 01 .009 06 


0.071 00 .161 07 0.009 32 
.072 91 2 .162 14 -009 59 
.074 90 ‘ . 163 22 -009 86 
.076 96 ‘ 164 31 .010 14 
.079 11 : .165 40 .010 42 


0.081 34 .166 50 0.010 72 
.083 66 ‘ . 167 61 -011 02 
.086 08 ; - 168 73 -O11 33 
.088 59 ‘ -169 85 -O11 64 
.091 21 ‘ .170 98 -O11 96 


.172 12 0.012 30 
.173 27 .012 64 
.174 42 (012 98 
.175 58 .013 34 
.176 75 .013 71 


.177 92 0.014 08 
.179 11 .014 46 
.180 30 .014 86 
.181°50 -015 26 
.182 71 .015 67 


.183 93 0.016 09 
.185 15 .016 52 
.186 38 -016 96 
. 187 62 .O17 41 
.188 87 .017 87 


.190 12 0.018 35 
.191 38 .018 83 
.192 65 .019 33 
.193 93 .019 83 
.195 21 -020 35 





s*/(% — xo)* 


0.196 51 
.197 81 
.199 12 
.200 43 
.201 75 


0.203 09 
.204 43 
.205 77 
.207 13 
.208 49 


0.209 86 
.211 24 
.212 62 
.214 01 
.215 41 


0.216 82 
.218 24 
.219 66 
.221 02 
.222 53 


0.223 98 
.225 43 
.226 89 
.228 36 
.229 84 


0.231 32 
.232 81 
.234 31 
.235 82 
.237 33 


.238 85 
.240 38 
.241 91 
.243 45 
.245 00 


0.246 56 
.248 12 
.249 69 
.251 27 
.252 85 
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TaBLE 1—(Continued) 
Auziliary estimation function, 6, for singly truncated samples 


0.020 88 
.021 42 
.021 98 
.022 54 
.023 12 


0.023 72 
.024 32 
.024 94 
.025 57 
.026 22 


0.026 88 
.027 55 
.028 25 
.028 95 
.029 67 


0.030 41 
.031 16 
.031 93 
.032 72 
.033 52 


0.034 33 
.035 17 
.036 02 
.036 89 
.037 78 


0.038 69 
.039 62 
.040 56 
.041 53 
.042 51 


0.043 52 
.044 54 
.045 59 
.046 65 
.047 74 


0.048 85 
.049 98 
.051 14 
.052 31 
.053 51 


s?/(% — x0)* 


0.254 44 
.256 04 
.257 65 
.259 26 
.260 88 


.262 50 
.264 14 
.265 78 
. 267 42 
.269 07 


0.270 73 
.272 40 
.274 08 
.275 74 
.277 43 


0.279 12 
.280 82 
.282 52 
.284 23 
.285 94 


0.287 66 
.289 39 
. .291 12 
.292 86 
.294 60 


0.296 35 
.298 11 
.299 70 
.301 63 
.303 40 


0.305 18 
.306 98 
.308 75 
.310 54 
.312 34 


0.314 14 
.315 95 
.317 76 
.319 57 
.3821 40 


0.054 73 
.055 98 
.057 25 
.058 55 
.059 87 


0.061 21 
.062 58 
.063 98 
-065 41 
.066 86 


0.068 33 
.069 84 
.071 37 
.072 93 
.074 52 


0.076 14 
.077 79 
.079 47 
-081 18 
.082 92 


0.084 69 
-086 49 
.088 33 
.090 20 
.092 10 


0.094 03 
.096 00 
.097 99 
-100 0 
.102 1 


.104 2 
.106 4 
-108 5 
.110 8 
.113 0 


.115 3 
.117 6 
.120 0 
.122 4 
.124 9 
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TaBLE 1—{Continued) 
Auziliary estimation function, 6, for singly truncated samples 


8*/(% a Xo)? 87/(Z = xo)? 


9.323 23 : 0.399 21 
.325 06 a -401 16 
.326 90 : -403 11 
.328 73 ; -405 07 
.330 57 d -407 02 


0.332 43 : 0.408 97 
.334 28 ; -410 90 
.336 13 : .412 88 
.338 00 ; .414 83 
.339 86 , .416 80 


0.341 73 : 0.418 76 
.343 61 : .420 72 
.345 48 ; -422 67 
.347 36 ‘ -424 63 
.349 25 ‘ .426 59 


0.351 13 : 0.428 53 
.353 02 ; .430 51 
.354 92 ; .432 47 
.356 82 : .434 43 
.358 72 : .436 39 


0.360 62 é 0.438 35 
.362 53 ; .440 31 
.364 43 é 442 27 
.366 35 d .444 23 
.368 26 ; .446 19 


0.370 18 0.448 15 
.372 10 ; .450 10 
.374 02 : .452 06 
.375 95 é .454 02 
.377 88 : .455 97 


0.379 81 ‘ 0.457 92 
.381 74 ‘ .459 88 
.383 67 . .461 83 
.385 61 i .463 78 
.387 55 ‘ .465 73 


0.389 47 ; 0.467 67 
.391 43 ; .469 62 
.393 37 F .471 57 
.395 32 , .473 51 
.397 27 : .475 45 
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TaBLE 1—(Continued) 
Auxiliary estimation function, 0, for singly truncated samples 


8?/(% — 2)" 8*/(Z — Xo) 


0.477 39 0.552 83 
479 32 
.481 26 
.483 20 
.485 13 


0.487 06 
-488 99 
-490 91 
-492 84 ‘ : 
.494 76 .569 02 


0.496 68 0.570 80 
.498 63 .572 63 
.500 51 : .574 34 
.502 42 .576 10 
.504 33 .577 86 


0.506 28 0.579 65 
.508 14 a .581 36 
.510 04 ; .583 11 
-511 93 .584 85 
.513 85 .586 58 


0.515 72 0.588 31 
.517 61 : .590 04 
.519 49 ; ..591 76 
.521 38 > .593 47 
.523 25 : .595 18 


0.525 13 0.596 89 
-526 97 ; .598 59 
.528 87 : .600 28 
-530. 74 : .601 97 
.532 60 .603 66 


0.534 46 0.605 34 
.536 31 .607 01 
.538 16 .608 68 
.540 01 .610 35 
.541 85 .612 O01 


0.543 69 0.613 66 
.545 53 d .615 31 
.547 36 3 .616 96 
.549 19 : -618 59 
.551 01 : .620 23 
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TasBLe 1—(Continued) 
Auziliary estimation function, 6, for singly truncated samples 


s/(& — ao} #2 — ae¥ 


0.621 86 : 0.777 12 
-623 48 5 -781 95 
.625 09 . -786 66 
.626 71 : -791 26 
-628 31 ‘ .795 74 


0.629 91 7 0.800 12 
-631 51 i .804 39 
.633 10 ‘ .808 55 
.634 68 ‘ -812 62 
.636 26 -816 58 


0.637 84 ; 0.820 44 
-639 40 , -824 21 
.640 97 , .827 88 
.642 52 : .831 47 
.644 08 3 .834 96 


4’ e@e® 


0.645 62 i 0.838 37 
.647 35 ‘ .841 69 
.648 70 . -844 93 
.650 23 : .848 09 
.651 75 ; .851 17 


BRERS RSBse 


eons 


0.653 27 ; 0.854 17 
.660 76 ; .857 10 
.668 14 , -859 95 
.675 36 ‘ .862 74 
.682 44 , .865 45 


OOO Ow 
SRREa 


0.689 38 k 0.868 
.696 18 é 871 
-702 84 : .873 
.709 36 R 876 
.715 74 : .878 


0.721 98 , 0.880 
728 08 
. 734 05 
.739 88 
.745 58 


0.751 16 
.756 60 
.761 91 
.767 10 
772 17 
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FIGurRE 1—Estimation curve for singly truncated samples 
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Although estimators (13) have been derived for samples that are singly 
truncated on the left, they are equally &’pplicable when samples are singly irun- 
cated on the right, as a consequence of the symmetry of the normal probability 
density function (1). In both cases @ > 0, and as shown in the sketch below, when 
truncation is on the left, 


(€— 2%) >0, 
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and 
a<Z whereas when 


truncation is on the right 


(z a Xo) < 0, 
and 


go sh 


Truncation on the Left at xo Truncation on the Right at xo’ 


On the lower scales, ¢ is the standardized normal deviate, t = (x — u)/o. Thus, when x 
is normal (yz, a), then ¢ is normal (0, 1), and if zr» = — zo’ , it follows that — = —?’. 


3. SmInGLY CENSORED SAMPLES 


We consider two types of censored samples. A Type I Censored Sample is 
one in which the terminus or point of censoring is fixed, while a Type II Censored 
Sample is one in which the number of censored observations is fixed. Within 
each of these categories, we may have censoring either on the right or left. Here 
as in the case of truncated samples, the symmetry of f() makes it unnecessary 
to consider both left and right censoring in detail, and again our derivations are 
confined to left restricted samples. 


Type I Singly Censored Samples 


In this category, we consider samples consisting of a total of N observations 
subject to the restriction that full measurement (i.e. unrestricted observation) 
of the random variable z is possible if and only if 


(a) x > 2 , in which case censoring is on the left, 


(b) x < x, in which case censoring is on the right, 


where 2p is a known fixed terminus. Let n designate the number of fully measured 
observations and n, the number of censored observations for which it is known 
only that « < 2» (or x > 2 , in the case of censoring on the right). Since 2» 
and N are fixed, both n and n, are random variables subject to the condition 
that n, +n = N. 

The likelihood function for a sample of this type is 


P = [N/m in!][F@I"(o V2) exp |- z (x; — 1/20 |, (15) 


where ¢ and F(£) are given by (2). 
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In this case the maximum likelihood estimating equations are 


Yo — wp = at, 
E—p=cY, 
s* + (¢ —.u)’ = o*[1 + EY], 
where 
Y(h, &) = [h/1 — A)JZ(-, with h=n,/N. (17) 


The first equation of (16) comes from (2) and is identical with the first equation 
of (4) for the truncated case. The last two equations of (16) result from taking 
logarithms of (15), differentiating with respect to u and o” in turn, and equating 
to zero. Here as in the truncated case, Z and s’ are the sample mean and variance 
respectively. 

Estimating equations (16) which apply in the censored case differ from equa- 
tions (4) which apply in the truncated case only in that Z(&) appearing in (4) 
has been replaced in (16) by Y(h, €) as defined by (17). Procedures analogous 
to those employed in the truncated case enable us to replace the system of 
equations (16) with the equivalent system 


o=s' + XE — XH)’, 
w= €— NE — %), 
fli-— Y(Y -—O)|Y — * = #/(@ — 2)’, 
where 
Ah, &) = Y(h, &)/[V(h, &) — €]. (19) 


Let £ designate the solution of the third equation of (18), let } = ACh, &) and 
the desired estimators become 


= + NE re to)’ 
p= & — ME — x). 


(20) 


As computational aids in this case, tables and graphs of A as a function of h 
and [1 — Y(Y — #]/(Y — &)? have been prepared. Since & is the solution of 
the third equation of (18), we determine \ directly for any given sample as 
that value of \ which corresponds to the sample statistics h and s’/( — 2)’. 
As in the truncated case, only one table is required. 

With h, x , Z, s’ and therefore s’/(# — 2,)” available from the sample data, 
it is necessary only that i be read from Table 2 or from the graphs of Figures 
2 or 3 as required, and that ¢° and f be calculated using estimators of (20). In 
Figure 2, \ is graphed for h = 0(.01).27, while in Figure 3, it is graphed for 
h = 0(.05).75. In both Figures 2 and 3, the range of s’/(Z — 2)” is (0, 1.3). 
In Table 2, \ is given to 4D for h = .01(.01).05(.05).50 and for s’/( — 20)” 
0(.05)1.00. 

As in the truncated case, £ can when required, be obtained from (14). Estima- 
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Figure 2—Estimation curves for singly censored samples h = 0(.01).27 


tors (20) are equally applicable to both left and right censored samples for the 
same reasons that estimators (13) apply to both left and right truncated samples. 


Type II Singly Censored Samples 


In this category, we can consider samples consisting of N observations of a 
random variable with probability density function (1) such that 


(a) the smallest N — n observations are counted but not otherwise 
measured (in which case censoring is on the left), 


(b) the largest N — n observations are counted but not otherwise 
measured (in which case censoring is on the right). 


Let x, designate the smallest (or largest) completely measured observation, 
and the sample thus consists of n completely measured observations each of 
which is equal to or greater than x, (or equal to or less than z,) plus N — 7 
unmeasured observations about which it is known only that « < x, (or z > 2,). 
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Figure 3—Estimation curves for singly censored samples h = 0(.05).75 


Estimators for this case turn out to be identical with those of (2) for Type I 
Singly Censored Samples, when we let 


Ln = Xo, (21) 
N-n=n,. 


Although there are no essential differences between estimators for Type I 
Singly Censored Samples and those for Type II Singly Censored Samples, 
variances of these estimators differ in the two cases as is shown in Section 4. 


4. SAMPLING ERRORS OF ESTIMATES 


The asymptotic variance-covariance matrix of (f, ¢) is obtained by inverting 
the matrix whose elements are negatives of expected values of the second order 
derivatives of logarithms of the likelihood functions. Accordingly we obtain 
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V(a) ~ [6°/E(n) N1@22/(@11422 — $%2)], 
V(é) ~ [8/E@)]G11/Gr1622 — $12), (22) 
Cov (a, é) ~ [6/E(n) [—412/@11G22 — 12)1, 


where E(n) is the expected value of the number of completely measured observa- 
tions, and gi: , Yi2 , G22 are respectively —[o”/E(n)|E(@°L/dp’), —[o°/E(n)] 
E(@?L/du 8c), and —[o’/E(n)|E(é°L/dc’). To simplify the notation, L has 
been written for In P, and ¢,; for ¢,;(£) or ¢;;(h, £) as applicable. 

In truncated samples and in censored samples of Type II, n is fixed, therefore 
E(n) = n. In samples that are Type I singly censored on the left, E(n) = 
N{1 — F(é]. The ¢,; for singly truncated and for singly censored samples of 
Types I and II are 


Truncated Samples Type I Censored Samples 
gut) = 1 — Z®@[ZE) — é], gut) = 1 + ZE)[Z(—é) + €], 
git) = Z(E){1 — EZ) — €]}, — pral€) = ZE){1 + E[Z(—£) + €]}, 
Gort) = 2 + Epil), goolt) = 2 + Evil), 


Type II Censored Samples 
gulh,&) = 1+ Y(h, [Z(—é) + €], 
gis(h, t) = Y(h, £){1 + &[Z(—é) + él}, 
Paalh, t) = 2 + Epirlh, £). 


It is to be noted that as N — ~, theg,; for type II censored samples approach 
the y,; for type I censored samples. Likewise as N > ©, then [1 — F(£)] > 
limy_. (n/N). Therefore, limiting values of estimate variances and covariance 
for both types of censored samples in this sense become equal. 

The variance and covariance of estimators based on samples that are restricted 
on the left may be calculated by substituting appropriate values of ¢,;(£) from 
(23) into (22), where £ = (x) — a)/é or £ = (x, — @)/¢ as applicable. For samples 
that are restricted on the right, calculations are the same except that ¢;;(— £) 
from (23) are substituted into (22). 


5. ILLUSTRATIVE EXAMPLES 


To illustrate the practical application of estimators derived in the preceding 
sections of this paper, one example of each of the six sample types considered 
here has been selected. 

Example 1. Left truncated. To insure meeting a lower specification limit of 
0.1215 in. on the thickness of a certain insulating washer, all production of this 
component is sorted through go, no-go gages, and all of thickness less than 
this value are discarded. For a random sample of 100 washers selected from the 
screened (i.e. the retained) production, it is found that ¢ = 0.124624 and s° = 
2.1106 X 10°°. Since n = 100 and z, = 0.1215, then (@ — x) = 0.003124 and 
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s’/(& — 2)” = 0.21627. By linear interpolation in Table 1, we obtain § = 0.03012. 
(Even without Table 1, this value might have been read from Figure 1 to three 
decimal places as 0.030.) Under the assumption that x is normally distributed, 
we employ estimators (13) and calculate 


6? = 2.1106 XK 10°° + 0.03012(.003124)” = 2.405 x 107°, 


and 
é = 0.00155, 


fp = 0.124624 — 0.03012(.003124) = 0.1245. 
From (14) 


£ = (0.1215 — 0.1245)/.00155 = —1.94. 


In determining the asymptotic variances and the covariance of f and 4, 
Z(—1.94) = 0.062399 is calculated from the defining relation of (5) with the 
aid of ordinary tables of normal curve areas and ordinates. This value might 
have been obtained from “The normal probability function: Tables of certain 
area-ordinate ratios and their reciprocals,” published as an editorial in Bio- 
metrika, Vol. (42), (1955), pp. 217-22. Using the truncated sample formulas of 
(23), we calculate ¢,,(—1.94) = 0.8751, ¢:2(—1.94) = 0.3048, and ¢..(—1.94) = 
1.4087. Using these values with E(n) = n = 100, and with ¢’ as calculated above, 
we employ (22) to calculate V(a) ~ 2.98 X 107°, V(a) ~ 1.85 X 107°, and 
Cov (a, ¢) ~ —0.65 X 107°. It then follows that 


og = VV) ~17X10", & = VV) ~14 xX 10%, 


and 
pi.e = Cov (A, 6)/V V(a)V(4) ~ —0.28. 


In the case of truncated samples and type I censored samples, these calcula- 
tions might be somewhat simplified with the aid of tables of elements of the 
variance-covariance matrices given by Hald [8]. Hald’s tables are also available 
in his “Statistical tables and formulas”, published by John Wiley and Sons 
(1952). Similar tables for type I censored samples were given earlier by Stevens 
[14]. Gupta [7] tabled corresponding matrix elements for type II censored 
samples, while the author and Woodward [2] tabled the matrix element necessary 
for calculating V(¢) in the truncated case. 

Example 2. Right truncated. To insure meeting a maximum weight specification 
of 12.00 ounces on a certain radio component for an aircraft installation, all 
production of this component is weighed and those units which exceed the 
specified maximum weight are discarded. For a random sample of 50 units 
selected from the screened (accepted) production, = 9.35 ounces, s” = 1.1264, 
n= 50, and 1.= 12.00. Accordingly (¢ — x,) = —2.65, and s°/(@ — xo)” = 0.16040. 
By interpolation from Table 1, 6 = 0.00916 (a value which might have been 


read to three decimals from Figure 1 as 0.009) and estimators (13) enable us to 
calculate 
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6” = 1.1264 + 0.00916(—2.65)? = 1.1907, 


é = 1.091 ounces, ; 
f = 9.35 — 0.00916(—2.65) = 9.38 ounces. 


From (14). £ = (12.00 — 9.38)/1.091 = 2.200. 

Since this sample is truncated on the right, we need the ¢;;(—£) in order to 
determine asymptotic variances and the covariance of these estimates. Accord- 
ingly, from (5) we calculate Z(— 2.200) = 0.035975, and from (23) ¢,,(—2.200) = 
0.9196, ¢12(—2.200) = 0.2129, ¢..(—2.200) = 1.5315. As in example 1, we 
substitute these values with E(n) = n = 50 and with é” as determined above 
into (22) to calculate V(f@) ~ 0.027, V(¢) ~ 0.016, and Cov (f, ¢) ~ —0.004. 

Example 3. Left Censored Type I. A certain breaking strength test is performed 
by applying an initial (minimum) stress and then gradually increasing it until 
the test specimen fails. To save time, the initial stress x) , may be established 
high enough that an occasional specimen fails under this minimum stress. For 
such censored readings, it is known only that the breaking strength is less than 
or at most equal to x, . Individual readings are obtained on all specimens for 
which x > x, . The resulting samples are thus singly censored on the left at x, . 
There is, however, one minor difference between samples involved here and those 
discussed in Section 3. Here the terminal is included in the interval where 
censoring occurs. In Section 3, the terminal was included in the interval of 
measurement. The maximum likelihood estimating equations turn out to be 
identical in the two cases so the difference is not important. For a hank strength 
test of the above type performed on a woolen yarn with an initial stress of 2, = 
70.00 lbs being applied, = 80.654, s’ = 31.73488, n = 50, and n, = 3. Accord- 
ingly ( — x») = 10.654, s’/(—2»)” = 0.27958, and h = 3/(50 + 3) = 0.056604. 
Two way interpolation from Table 2, gives } = .071, and from estimators (20) 


6° = 31.73488 + .071(10.654)” = 39.72393, 
and 


é = 6.31, 
p = 80.654 — .071(10.654) = 79.90, 
— = (70.00 — 79.90)/6.31 = —1.57. 


For a more accurate determination of these estimates, we calculate additional 
values of Y, \, and related functions from the defining relations given in Sections 
2 and 3 with the aid of normal curve tables of ordinates and areas or from the 
Biometrika editorial tables (loc. cit.). We then interpolate as summarized 
below to obtain i to additional significant digits. 


s’°/( — 4%)” d g 


.28347 .07109 — 1.560 
27958 07096 —1.569 
27921 07095 —1.570 
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With \ thus more accurately determined as .07096, we again employ (20) to 
calculate final estimates as 







6? = 39.79365, é¢ = 6.308, 
p = 79.90, and £ = —1.569. 


To determine asymptotic variances and the covariance of these estimates, 
we compute Z(—1.569) = 0.12372 and Z(1.569) = 1.99755. Using the type I 
censored sample formulas from (23), we compute ¢,,(—1.569) = 1.0532, 
¢i2(— 1.569) = 0.04053 and ¢..(—1.569) = 1.93641. Using these values with 
E(n) = 53[1 — F(—1.569)] = 49.91, and with é* as computed above, we employ 
(22) to calculated variances and the covariance. Finally, o; = V V(a) ~ 0.87, 
o; = VV(é) ~ 0.64 and p;., ~ —0.03. 

Example 4. Right Censored Type I. A reaction time test is terminated at the 
end of ten hours in order to eliminate the effects of certain contaminants which 
are troublesome when the test is continued over a longer period. For a specific 
sample of this type, x) = 10, n = 62, n, = 38, = 8.75, s” = 1.1043, ( — x) = 
—1.25, s°/( — 2)” = 0.70675, and h = .38. Two way linear interpolation in 
Table 2 gives } = 0.71. Accordingly, using estimators (20) we calculate 













8.75 — .71(—1.25) = 9.64, 
1.1043 + .71(1.5625) = 2.2137, 
1.49, and £ = 0.244. 


Bp 
3? 


























é 


Since this sample is censored on the right, we need the ¢:i(—8. in order to 
determine the variances and the covariance of @ and ¢. Accordingly, we calculate 
values of Z(—0.244) and Z(0.244) as defined by (5). From the type I censored 
formulas of (23) we evaluate ¢,,(—0.244), ¢2(—0.244), and 22.(—0.244). 
With E(n) = 100[1 — F(—0.244)] = 100F (0.244), and ¢’ as calculated above, 
we employ (22) to calculate V(a) ~ 0.070, V(é) ~ 0.301, and Cov (f, ¢) ~ 
—0.132. 

Example 5. Left Censored Type II. In order to estimate the mean and standard 
deviation of the distribution of weights of a certain manufactured component, 
a random sample of size 60 is selected. By means of a rapid weight-sorting 
procedure, the 5 lightest specimens are identified and no further measurements 
are made on them. Thus, it is known only that each of these 5 censored specimens 
weighs less than z, , which is the weight of the lightest of the 55 remaining 
specimens. Each of these remaining specimens is weighed and for the sample 
selected € = 16.35 ounces, s° = 8.643, n, = 5,n = 55, and x, = 9.27. Accordingly, 
h = 08333, (€ — x,) = 7.08, and s’/(é — z,)” = 0.17242. Two way linear 
interpolation from Table 2 gives } = .10 and using estimators (20) we calculate 


pf = 16.35 — 0.10(7.08) = 15.64, 
é° = 8.643 + 0.10(50.1264) = 13.656, 
é¢ = 3.69, and £& = (9.27 — 15.64)/3.69 = —1.73. 
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To determine asymptotic variances and the covariance of sample estimates 
in this case, we use (5) to calculate Z(1.73) = 2.13637, and (17) to calculate 
Y (0.08333, —1.73) = [5/55]Z(1.73) = 0.194215. Using these values in the type 
II censored formulas of (23), we calculate ¢,,(0.08333, —1.73) = 1.07892, 
%12(0.08333, —1.73) = 0.02344, and ¢,2(0.08333, —1.73) = 1.95945. Since 
this is a type II censored sample, E(n) = n = 55. Using these values and that 
of é° as determined above, we employ (22) to calculate V(a) ~ 0.230, V(¢) ~ 
0.127, and Cov (f, ¢) ~ —0.003. 

Example 6. Right Censored Type II. A sample of N = 300 electric light bulbs 
were life tested until nm = 119 had burned out with the result that = 1304.832 
hrs., s” = 12128.250, and x, = 1450.000 hrs. Accordingly s*/(@ — 2,)’ = 0.575515, 
nm, = 300 — 119 = 181, and Ah = 181/300 = 0.6033. Visual interpolation from 
the graph of Figure 3, gives } = 1.36, and using estimators (20) we now calculate 


A = 1304.832 — 1.36(1304.832 — 1450.000) = 1502 hrs., 
6? = 12128.250 + 1.36(1304.832 — 1450.000)” = 40789, 
and 
202 hrs. 
From (14) 
£ = (1450 — 1502)/202 = —0.257. 


This example was originally given by Gupta [7], and to the number of significant 
digits given, the above estimates are in agreement with those which he calculated. 
A more accurate determination of { and correspondingly more accurate deter- 
minations of @ and ¢ are possible by calculating additional values of A}, Y and 
related functions directly from tables of the normal curve areas and ordinates 
or from the Biometrika editorial tables (loc. cit.) as in example 3, and then 
interpolating. Following is a summary of such calculations for the above data. 


s°/(é =a 2)” —— A 


0.575304 0.25690 1.35712 

0.675515 0.25693 1.85719 

0.576081 0.257: 1.35735 
With } = 1.35719 as determined above, a recalculation using estimators (20) 
gives more accurate values as @ = 1501.853, ¢ = 201.815, and of course £ = 
— 0.25693. 

This is a right censored type II sample, and in order to determine variances 
and the covariance of the sample estimates, we must evaluate the ¢,;(h, —8; 
that is, ¢1:(h, 0.257), g:2(h, 0.257) and g22(h, 0.257), where h = 0.6033. Calcu- 
lating these values in the same manner as for example 5, with E(n) = n = 119, 
and with 6” as determined above, we substitute in (22) and subsequently calculate 
o, = VV(a) ~ 16.6, 0, = V V(6) ~ 14.39, and p;,, ~ —0.57. The rather high 
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correlation between estimates here, reflects the high degree of censoring in this 
example. 

The assistance of Mr. Walt G. Horstman, who performed the calculations 
necessary for the compilation of Tables 1 and 2, and who rendered material 
aid in the preparation of the charts of Figures 1, 2, and 3 is gratefully acknowl- 
edged. 
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Control Chart Tests Based on 


Geometric Moving Averages 


S. W. Roserts 
Bell Telephone Laboratories 


A geometrical moving average gives the most recent observation the greatest 
weight, and all previous observations weights decreasing in geometric progression 
from the most recent back to the first. A graphical procedure for generating geometric 
moving averages is described in which the most recent observation is assigned a 
weight r. The properties of control chart tests based on geometric moving averages 
are compared to tests based on ordinary moving averages. 


1. INTRODUCTION 


An X control chart presents a graphical summary of information collected for 
the purpose of providing continuing information on the fv: :tioning of a process. 
This information is useful in indicating the existence, and in some instances the 
causes, of conditions requiring corrective action. A central line is drawn on the 
control chart to represent the level at which the process operates in the absence 
of trouble, and limit lines are placed at prescribed distances above and below 
this central line. The presence of trouble is indicated if the point X; , repre- 
senting the average of several measurements taken at time j (j = 1, 2, 3, ---), 
falls outside either of the limit lines. (Symbols such as X; are used to denote both 
a point and its scale value). In practice, corrective action is sometimes taken when 
arun of points occurs that is highly improbable in the absence of trouble, though 
no single point falls outside the limit lines; for example, a run of seven consecu- 
tive points on one side of the central line is sometimes taken as an indication 
that trouble is present. 

In a recent paper concerned primarily with the properties of run tests in 
detecting shifts in a process average, (1) the writer indicated that tests based on 
moving averages should be more efficient than run tests in making use of the 
information supplied by control chart points. A graphical procedure was outlined 
for generating points representing moving averages of 2, 4, 8 or 2° consecutive 
points. The procedure was based on the fact that the point half-way between 
two points on their connecting straight line represents their average value. For 
example, points representing the average of eight X’s are obtained from points 
representing two sets of four X’s, which in turn are obtained from points repre- 
senting pairs of X’s. The necessity of generating the auxiliary points representing 
four and two X’s is a disadvantage; such points could be incorporated into the 
testing procedure by drawing limit lines for them but it is questionable whether 
the improvement in statistical properties would justify the added complexity. 
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J. W. Tukey amit to the writer an alternative graphical procedure that 
requires no auxiliary points for generating a type of moving average in which the 
entire history of X’s is assigned weights, with weights decreasing as a geometric 
progression from the most recent point back to the first. Such moving averages 
are called geometric moving averages in this paper, while the term moving average 
without the «jective “geometric” refers to the ordinary moving average in 
which k consecutive points are assigned equal weights of 1/k. In computing a 
geometric moving average the most recent point is assigned a weight of r, 
0 <r < 1, and every other point is assigned a fraction (1 — r) of the weight of 
its immediate successor. 

This paper outlines a procedure for generating geometric moving averages, and 
compares the properties of control chart tests based on such averages with tests 
based on ordinary moving averages. 


2. GENERATING GEOMETRIC Movine AVERAGES 
Let Z; denote the geometric moving average at time j (j = 0, 1, 2, 3, ---). 

It is convenient to set Z) = po , where yp is the ordinate of the central line of 
control chart. Thus, 

Zo = bo ’ 

Z; (1 a r)Z;-1 + rX; ’ j > 0. 
The point Z; is located a fraction r of the distance from Z;_, to X; on their con- 
neciing straight line. 


To exhibit the geometric progression of weights assigned to the X’s, (1) can 
be written as 


(1) 


Z, = (1 —1)'uo + (1 — 1)" RX, +71 — 1) Xs 
+--- +1 —9F,,.492,, § 20. 
Rearranging terms, this becomes 


(Z; — mo) = r(l — r)i-"(X, — to) + (1 — r)'*(X, — Ho) 


(2) 


(3) 
+ ee) $l — (Xj. — wo) +17(X; — wo), 7§ > 09, 


and clearly the expected value of Z; is 49 when the expected values of the X’s 
are all equal to uy . Then, if the X’s are independent and have a common standard 
deviation og , the standard deviation of Z; is 


oz, = ds —[1-(-n*ler, (4) 


which increases to its limiting value 
oz = Vr/(2 — 1) og 
as j increases. 
To generate the Z’s, proceed as follows: 


1. Mark point X,.at abscissa 1. 
2. Mark Z, = uo at abscissa 1 — 1/r. 
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3. Draw a straight line connecting Z,) with X, ; the point on this line with 
abscissa 1 — (1 — r)/ris Z, . 


The next point, Z, , will be at abscissa 2 — (1 — r)/r on the straight line con- 
necting Z, with X, , which is plotted at abscissa 2. Both the X’s and the Z’s 
progress one abscissa unit at a time; Z; is always (1 — r)/r abscissa units to the 
left of X; . 

Figure 1 illustrates the procedure for the case where r = 3, which is a value of 
r that may have wide appeal in practice. For one thing, with r = 2 the Z’s 
fall on vertical lines drawn half-way between vertical lines of the X’s, thus 
minimizing confusion between the two types of points. Also, with r = 2, ¢z = ox/2 
so that, as shown in Figure 1, 3c limit lines for the Z’s are half as far from the 
central line as 3 limit lines for the X’s. 
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Ficurs 1—Generating geometric moving averages. 
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Limit lines for the Z’s should be based on az ;i.e., place limit lines at up + 307 . 
In general it wouldn’t pay to widen the limits with increasing j using equation (4). 

The sequence of Z’s may be started at any time j simply be setting Z;_1 = uo 
and following the procedure outlined above. Ai: error ¢ in Z;_, manifests itself 
as an error e(1 — r)* in Z;-44: . 


3. Process MopEL 


One way of comparing control chart procedures is to assume that the end 
result is a sequence of decisions on whether or not to take action. An evaluation 
of the effectiveness of various tests must center attention on the relative fre- 
quencies of two types of decision errors, with their attendant economic conse- 
quences. In discussing the two types of decision errors it is convenient to adopt 
the terminology of hypothesis testing and define a null hypothesis; then erroneous 
rejections and acceptances of the null hypothesis are called, respectively, Type 1 
and Type 2 errors. 

It is assumed that at any time there is a process average » = wo + 5, where 
Mo is the nominal process average and 6, which is measured in the same units 
AS Uy , is the difference between the current process average and its nominal value. 
In the absence of trouble 6 = 0. 

The testing procedure consists of periodically testing the null hypothesis that 
= uo, Or, equivalently, that 5 = 0. A rejection of the null hypothesis when it is 
“true” is a Type 1 error; such an error results in wasted effort of looking for 
nonexistent trouble. An acceptance of the null hypothesis when it is not true is 
a Type 2 error; the cost of such an error must be expressed in terms of the rela- 
tively poor quality of output from time j to time j + 1, when the null hypothesis 
is next tested. 

The over-all average quality of the product depends on the way the process. 
average varies with time, as affected by occurrences of trouble and by attempts. 
to keep the process free of trouble. A simplifying concept is that the process 
average remains at its nominal level until trouble shifts it to a new level where 
it remains until the trouble is eliminated. 

Let the random variable y be the number of points following the occurrence. 
of trouble until the decision is made that trouble is present. A sequence of y — 1 
consecutive Type 2 errors is made. The distribution function of y depends not 
only on the amount 6 of the shift but also on the standard deviation of X. For 
this reason A = 6/cg is used, rather than 6, as a measure of the amount of the 
shift. 

For most practical purposes the distribution function of y is adequately 
summarized by its expected, or average, value E(y). Curves of E(y) versus A 
are presented in Figures 2-5 to describe the effectiveness of particular tests in 
detecting shifts in a process average. The test notation is defined in Section 4. 

Since X; is the average of n measurer.ients taken at time j, the assumption is 
made that its distribution is normal with average » = uo + 6 and standard 
deviation cg = o’/~/n, where o’ is the standard deviation of the individual 
measurements and is assumed to be known from past experience and not subject 
to change. In considering a curve of E(y) versus A for a particular application,. 
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FicurE 2—Moving average tests, M; . 


note that A = 6-~Vn/o’; thus, while the curves are independent of sample size n, 
the value of A applicable in a particular situation depends on n. 

For a particular value of A, the corresponding value of E(y) — 1 is the expected 
number of consecutive Type 2 errors made following a shift of the magnitude 
associated with A. The value of Z(y) at A = 0 is a measure of the frequency of 
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Figure 3—Moving average tests, Mi... 


Type 1 errors. For example, the top (for central values of A) curves in Figures 
2-5 pertain to the standard control chart test with 3c limits, and on these top 
curves the value of E(y) at A = 0 is roughly 370; this means that when the process 
average is at its nominal level, u , an average of 370 points can be expected 
between successive Type 1 errors when using the standard control chart test. 
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Fiaure 4—Geometric moving average tests, G, . 


4. Curves or E(y) Versus A 


Let M, denote the test that prescribes rejection of the null hypothesis at time 
j if the moving average of k consecutive X’s falls outside limit lines at up + 
3cx/~Wk. (If limits other than 3¢ limits are used for a statistic, the notation 
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Fiagure 5—Geometric moving average tests, Gi,r . 


will be supplemented to indicate the factor used—e.g., tests M, (2) employ 2c 
limits: uy + 20" Vk.) Let M, , denote a combination of this test with the standard 
control chart test, which prescribes rejection of the null hypothesis if X;, falls 
outside limit lines at uo + 3c . Similarly, let G, denote the test prescribing 
rejection of the null hypothesis if the geometric moving average with weight 
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parameter r falls outside limit lines at wo + 3Vr/(2 —nr ox , and let G,,, 
denote such a test in combination with the standard control chart test. Figures 2 
through 5 present curves of H(y) versus A for these tests for various values of 
the weight parameters. Note in each case that when the weight parameters 
r or k are given their limiting value, 1, the test becomes the standard control 
chart test; comparisons with the standard test are thereby facilitated. 

The curves are based largely on estimates of E(y) at A = 0(.5) 3.5 based on 
results of process simulations using 25,000 random normal numbers. Sets of 
eight numbers were used to represent the process prior to the occurrence of 
trouble; if the statistics computed with these eight numbers fell within their 
limits, an amount A was added to subsequent numbers and a count was made 
of the number y of such numbers until a statistic fell outside its limit. The average 
value of y for all such cycles (roughly 25,000/(E£(y) + 8) in number) serves as 
an estimate of E(y) for the particular test and A considered. Such estimates, 
other than at A = 0, fall on or very close to the curves shown. The estimates at 
A = 0 were adjusted to be consistent among themselves and with computed 
lower bounds. 

While the curves of Figures 2-5 afford rough comparisons of the statistical 
properties of various tests under the conditions assumed, more useful compari- 
sons could be made if the limit lines were set so that the tests had a common 
level of Type 1 errors. Such a comparison applied to the moving average tests 
of Figure 2 would show that each such test has its curve lower than the others 
in a particular interval of A, and these intervals approach A = 0 for increasing k; 
consequently, it would be well to neglect the history prior to the k most recent 
control chart points in applications where interest is in A greater than a particular 
value A, , which depends on k. 

With geometric moving average tests, the entire history of points is given some 
weight; by considering only the k most recent points some improvement could 
be effected in detecting changes greater than a critical value determined by k. 
However, the more remote history is given relatively little weight, and has 
correspondingly little influence on test properties. In fact, it appears that given 
any moving average test a geometric moving average test can be selected that 
has roughly equivalent properties; the requirement for rough equivalence appears 
to be that the o’s be equal, i.e., that r = 2/(k + 1). The choice between two 
tests having similar statistical properties must be based largely on their relative 
operational simplicity in the application under consideration. 

To illustrate the use of the curves, consider a resistor manufacturing process 
in which at the end of each hour four resistors are measured and the average 
of the measurements is plotted on an X control chart. Suppose experience 
indicates that 4» = 1000 ohms and o’ = 10 ohms. 

For test G, 2,5 the limit lines for the X’s are at 1000 + 3 X 10/4 ohms, 
or at 985 and 1015, and the limit lines for the Z’s are at 992.5 and 1007.5 (See 
Figure 1). Now suppose something happens to shift the process average to 1005 
ohms. How long will it be, on the average, before this shift is detected? 

Compute A = (1005 — 1000)/5 = 1, and look at Figure 5. There is no curve 
for r = 2, but an interpolation between values for r = } and r = } indicates that 
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Figure 6—Approximate differences between E(y) of various tests and E(y) of the standard 
control chart test; arithmetic scale. 


E(y) is about 14. Thus, an average of about 14 pairs of X’s and Z’s will be plotted 
before either (1) an X is greater than 1015 (or less than 985) or (2) Z is greater 
than 1007.5 (or less than 992.5). The top curve of Figure 5 indicates that an 
average of about 44X’s will be plotted before one of them is greater than 1015 
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(or. less than 985) when A = 1. Since the top curve pertains to the standard 
control chart test, it appears that test G,,2,5 will save an average of about 
44 — 14 = 30 hours in detecting a shift of five ohms in the process average in 
the example used here. 

On the other hand, Figure 5 shows that E(y) is about 250 at A = 0 for test 
G25 , a8 compared with about 370 for the standard control chart test (top 
curve). Using test G,2,s , then, when the process remains trouble-free an average 
of 250 hours will pass between false indications of trouble (Type 1 errors) as 
compared with 370 hours for the standard control chart test. 

Because the ordinate scales are logarithmic in Figures 2-5 it is difficult to 
appreciate the savings in E(y) for the various tests as compared to the standard 
control chart test. The top curve in Figure 6 shows roughly the savings afforded 
using either M,(2.95) or G,,,(2.95), with limit lines set so that these tests have 
approximately the same frequency of Type 1 errors as the standard control 
chart test. Clearly, the important savings occur for relatively small shifts, say 
| A | < 2, and the loss at large values of A is relatively small. 

The curve of Figure 6 called “Run Test” is for a test that prescribes rejection 
of the null hypothesis if (1) a single X falls outside limit lines at yp + 3.13¢¢ or 
(2) ten consecutive X’s fall on one side of the central line of the control chart. 
The curve represents roughly what can be achieved using multiple run tests— 
use of additional run criteria would raise the curve in some ranges of A and lower 
it in others, provided the limit lines were set to yield the same frequency of 
Type 1 errors as the standard control chart test so that the curve would go 
through the origin of Figure 6. 


5. CoNCLUSION 


Tests based on geometric moving averages compare most favorably with mul- 
tiple run tests and moving average tests with regard to simplicity and statistical 
properties. The standard control chart test is simpler than alternative tests 
and cannot be improved on in detecting relatively large changes. In fact the 
parameters that would be chosen for the other tests in such applications would 
reduce such tests to the standard control chart test. In general, as interest 
centers on early detection of smaller and smaller changes, the appropriate param- 
eter r of geometric moving average tests decreases from unity to smaller and 
smaller values, and the number k of consecutive points considered in moving 
average tests and multiple run tests increases. 

The complexity of geometric moving average tests relative to the standard 
control chart test makes its use difficult to justify in many cases. In other cases 
the added sensitivity to relatively small changes in the process is adequate 
justification. In still other cases the cost of the added complexity is more than 
covered by a reduction in sampling costs. For example, in some applications 
& geometric moving average test used on a chart where the averages or medians 
of three measurements are plotted will be more effective than the standard con- 
trol chart test based on the averages of four measurements. 

Sample size and frequency of sampling have been largely eliminated from the 
discussion, since interest has been primarily on making efficient use of the informa- 
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tion supplied by the control chart points. The preceding paragraph provides an 
example of the importance of considering these variables when selecting tests 
for particular applications. 

Geometric moving average tests would seem to be particularly well suited to 
applications where tests based on ordinary moving averages are now used, 
including cases where the sample size is unity. 

While this paper has been concerned specifically with X control charts, much 
of it could be rewritten word for word to apply to control charts of other sample 
statistics by simply replacing the symbol “X” and the corresponding term 
“average” with appropriate substitutes. Certainly the particular weight patterns 
applied to the history of points on an X control chart in obtaining ordinary 
moving averages and geometric moving averages can just as readily and meaning- 
fully be applied to the points on control charts of such useful sample statistics 
as the maximum, median, minimum and range. The curves of E(y) versus A 
are generally applicable if allowances are made for the non-normality of the 
distributions of the sample statistics—curves pertaining to tests which apply 
weights to several control chart points are affected very little by non-normality. 

It thus appears that the graphical procedure for applying geometric weight 
patterns may prove useful in a wide variety of applications. 
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This paper deals with the theory of a proposed method for the statistical study of 
measuring processes. The practical aspects of the method, including computational 
details, are discussed in a companion paper published in the ASTM Bulletin. In 
the present article a theoretical framework is proposed for the mathematical ex~- 
pression of the sources of variation in measuring methods and a suitable method 
of statistical analysis is described. Particular attention is given, both here and in 
the companion paper, to interlaboratory studies of test methods. An illustration 
based on data taken from the chemical literature is appended. 


INTRODUCTION 


This paper is concerned with the statistical study of methods of measurement. 
In order to avoid misunderstandings, we wish to distinguish clearly between 
two types of measurement. One type concerns the measurement of fundamental 
quantities, such as the acceleration of gravity, or the velocity of light. The 
second type of measurement relates to the determination of properties of 
materials, or of physical systems, such as the magnesium content of an alloy, 
or the impedance of an electronic piece of equipment under fixed conditions. 
The basic difference between these two categories is that in the first, the quantity 
measured is a unique constant, whereas in the second category, this quantity 
can assume a more or less extensive range of values. We will be concerned 
exclusively with measurements of the second type. 

No method of measurement is useful unless it satisfies certain minimum 
requirements of reproducibility. It is for this reason that considerable attention 
has been given, in recent years, to studies of test methods involving relatively 
large numbers of laboratories. These studies result from the often repeated 
observation that agreement between laboratories is, in general, much less 
satisfactory than can be achieved within a single laboratory, even in the case 
of well-established methods of measurement. Vast sums of money as well as 
large numbers of man-hours are being spent in studies of this nature. Naturally, 
attempts have been made to apply the principles of statistical design and 
statistical analysis to this subject. If, on the whole, the introduction of statistics 
in this field has not been so successful as might be expected, the reason may 
perhaps be found in the undesirability of attempting to force the subject of 
interlaboratory testing into the familiar frameworks of either factorial or 
hierarchical designs, without a critical examination of the true experimental 
situation underlying interlaboratory comparisons and without a clearcut state- 


251 





252 JOHN MANDEL 


ment of the practical goal of such studies in terms of the possible actions to be 
taken as a result of the analysis. 

The present paper is an attempt to reconsider the subject by placing it in a 
wider framework. A general model is proposed for the understanding of ex- 
perimental error in methods of measurement. The model is then somewhat 
particularized in order to make it possible to apply accepted methods of linear 
regression. The author believes, as a result of the examination of an appreciable 
number of interlaboratory test data, that the proposed linear model is satis- 
factory in most instances. (The practical aspects of the theory developed in 
this article, together with an application to a technological problem, are described 
in a companion paper [4].) It should be noted that the proposed model includes 
many of the more commonly considered models as special cases. 


Tue NATURE OF EXPERIMENTAL ERROR 


It is useful to look upon a method of measurement as a process involving 
three basic elements. The first element is the unknown quantity, Q, which is 
the object of the measurement. Generally, Q enters the measuring process in 
the form of a sample or specimen, or of a physical system. The second element 
is the numerical value, M, obtained by applying the measuring process, in 
accordance with a given set of instructions, to the specimen or the system. Of 
course, M is a function of Q, but even in the most precise methods of measure- 
ment there always exist factors other than Q, and related to environmental 
conditions, that have an effect on the measurement M. Let us denote the totality 
of these environmental factors by the letter C, keeping in mind that C repre- 
sents not a single variable, but rather a vector in a space with a large number 
of dimensions. The vector C is the third element involved in the measuring 
process. For any given value of the vector C, say C, , the measurement M 
can be considered to be a mathematical function of Q, M = f (Q | C,). It is 
assumed in this paper that this function is monotonic in Q over the entire range 
of validity of the method of measurement, since a non-monotonic function 
would be of little practical value. As C varies, there corresponds to each value 
of C a curve representing the relation between M and Q. The totality of these 
curves, for the entire range of C-vectors, constitutes the measuring process. 

For a well defined method, the factors constituting C vary over a relatively 
small range. The value taken by each of these environmental factors during a 
given measurement is to a certain extent governed by chance. Nevertheless, 
some of these factors are systematic in that they depend on identifiable entities, 
such as a given laboratory or instrument, a given day, a given operator, etc. 

Measurements corresponding to the same value of Q will not necessarily 
agree with each other, as they may lie on different M, Q curves. The variability 
thus produced is known as “experimental error’. It is seen that, depending 
on the extent to which the environmental conditions are controlled, experimental 
error may be large or small, and is partly systematic and partly random. 

When all identifiable components of C are “controlled”, i.e., when they are 
given fixed values to the maximum feasible extent, the remaining components 
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of C may be considered as random variables. Their variability generates a set of 
curves, which we may call a ““bundle’’, such that it is a matter of chance, which 
curve of the bundle prevails for any given measurement. In this manner there 
corresponds to each value of Q, a statistical population of M values, with a 
particular standard deviation, related to the “width” of the bundle for. this 
Q-value. This standard deviation may be a function of Q as well as of other 
factors such as, for example, the nature of the material or system measured. 
In practice it is generally found that the dependence of the standard deviation 
on Q sufficiently overshadows its dependence on all other factors to result in a 
recognizable functional relationship. The standard deviation thus described 
does not measure the full variability of the measuring process, since it does not 
include the variability due to the systematic factors for which control has been 
exercised. 


INTERFERING PROPERTIES 


So far, it has been tacitly assumed that for a given value of the vector C, a 
single curve expresses the relation between the measurement M and the magni- 
tude of the measured property Q. However, there exists probably no measuring 
process that responds to a single property of a portion of material or, more 
generally, of a physical system, at the exclusion of all its other properties. For 
example, in chemical analysis, one deals with the well-known problem of “inter- 
fering substances”. In physical measurements the situation is very much the 
same. For example, in the determination of the electrical conductivity of a 
material, the measurement is disturbed by the secondary phenomenon of surface 
conductivity, the variation of which, from one material to another, is not 
necessarily related to their volume conductivities. Of course, it may be possible 
in situations of this type, to make appropriate corrections for these disturbances. 
However, there are many situations, particularly in the field of technological 
measurements, where it is either impossible, or impractical to correct for all 
known disturbing factors, not to mention the influence of unknown interfering 
properties. Suppose now that an experiment is performed in which a number 
of specimens or systems covering the pertinent range of Q values, are subjected 
to the measuring process, under controlled test conditions. It then follows from 
the preceding observations that, apart from the experimental error caused by 
random variations in the C vector, an additional disturbance may be noticed, 
caused by the fact that the various samples vary not only in the property Q, 
but also in the interfering properties affecting the measured values. For clarity 
of presentation, we will refer to this type of disturbance as the I-factor, where 
I stands for “interfering properties’. More will be said about this factor, but 
it is worth mentioning, at this point, that some I-factors are such that their 
effect is virtually unrelated to the environmental conditions (value of C-vector), 
whereas others are highly sensitive to these conditions. Thus, in the example 
of conductivity measurements cited above, any environmental factor affecting 
the surface-condition of the test specimen is likely to have a far more noticeable 
effect on the surface-conductivity than on the volume-conductivity. Stated in 
statistical terminology, though perhaps somewhat loosely, this means that the 





254 JOHN MANDEL 
I and C factors may well have an appreciable interaction with respect to the 
measurement M. 

Inclusion of the I-factors leads to a more general representation of the measur- 
ing process, through the relation: 


(1) M = fQ, IC) 


This relation expresses the proposition that the measurement M, for any 
given value of the vector C (representing environmental conditions) is a func- 
tion of the property to be measured, Q, and of the interfering properties, I. 


THe EXPERIMENTAL SCHEME 


Theoretically, the study of a method of measurement would consist in the 
complete elucidation of equation (1). From a practical viewpoint, however, 
this is neither feasible nor desirable. In the first place, the quantity Q is seldom 
exactly known for any given specimen or system. Furthermore, the factors I 
and C encompass such a wide range that a selection must necessarily be made. 

While situations exist in which a series of specimens or systems can be pre- 
pared for which the value of Q is exactly known, such as for “synthetic samples’’ 
in analytical chemistry, we will not require this to be the case. It will be suf- 
ficient to be able to set up a series of specimens or systems, such that the values 
of Q corresponding to them cover the range over which the method is to be 
studied. 

In regard to the vector C, it will be assumed that the process of measurement 
has already been studied and specified to such an extent that the variables of 
real importance have been assigned optimum values insofar as this is practically 
possible. For example, if the process is a method of chemical analysis, such 
factors as size of sample, amount of reagents, duration of the various reactions, 
etc., insofar as they affect the method, are supposed to have been assigned 
reasonable values and these values are assumed to have been incorporated in 
the description of the method. Experience shows that despite such precautions, 
agreement of measurements made in different laboratories, or by different 
operators, or at different times, or on different instruments of the same type, 
often falls considerably short of expectation. Therefore, we will here be con- 
cerned with those C-factors that are associated with the variables just mentioned, 
or of a similar type. For convenience, our main classification factor will be 
referred to as “laboratories”, but this term should be taken in a broad sense. 
For example, in many cases, a “laboratory” may actually refer to a series of 
measurements made in a single session by a particular operator, using a par- 
ticular instrument, in a particular laboratory. 

Finally, with regard to the I-factors, a decision must usually be made at the 
start concerning the diversity of types of materials or systems to be covered 
by the study. For example, in a study of the determination of the smoothness 
of paper, a decision must be made regarding the inclusion of various types of 
paper, such as mimeobond, newsprint, paperboard, etc. For a thorough study, 
it is undoubtedly advantageous to conduct a separate study for each type. 
Against this must be weighed the scope of validity of the study as well as economy 
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Basic Design 


Material 


Laboratory 1 2 :-- j °° 
1 . 
2 


Yijk 


contains 
z,jcell |: n replicates 


factors. It should, however, always be borne in mind that the wider the scope 
of the study in regard to types of materials, the larger will be the disturbance 
caused by I-factors. This is, of course, especially true in chemical analyses 
where materials of different types contain different kinds of interfering sub- 
stances. 

Under the conditions and limitations that have just been discussed, a study 
of interlaboratory scope can be represented by the scheme shown in Table I. 
A selection of materials, b in number, and covering the range of interest, are 
measured in a “laboratories’’. Each laboratory makes n replicate measurements 
on each material. Let y;;, represent the kth measurement on material j by 
laboratory 7. The n replicates within any particular cell are considered as a ran- 
dom sample from a theoretically infinite population of measurements within that 
cell. On the other hand, the “laboratories” will not be considered as a random 
sample from a larger statistical population, but rather as a “fixed” variable 
(see Eisenhart [i]). Thus, inferences concerned with variability among labor- 
atories are confined, at least theoretically, to the participating laboratories. 
With regard to “materials” the set of Q-values corresponding to the b materials 
constitutes, of course, a “‘fixed’’ variable. But, as will be seen later, each material 
is considered as a random selection from a population of materials with the 
same @-value. 


THE STABILIZATION OF WITHIN-LABORATORY VARIANCE 


The plan in Table I lends itself to the evaluation of within-laboratory vari- 
ability by computing the standard deviation among the n replicates within each 
cell. Unless n is very large, an unlikely situation in practice, these standard 
deviations will show a good deal of variability. The question arises whether the 
“true” within-cell standard deviations (i.e., the “population” standard devia- 
tions rather than their sample estimates) are also different for the various 
cells. Experience shows that in the case of reasonably well established methods 





256 JOHN MANDEL 


of measurement, the within-cell variability is essentially the same for most 
laboratories. On the other hand, it often varies considerably from material to 
material. This variability may be due to two causes. The first is inherent in the 
design of the study, which demands that the materials cover an appreciable 
range in the magnitude of the measured quantity, Q. Since the standard deviation 
may well be a function of Q, one often detects a systematic change in within-cell 
standard deviations, as one proceeds from materials for which the measured 
quantity is small to those for which it is large. The second cause of variability 
of the standard deviation among materials is due to differences in nature 
between the various materials, resulting in differences in the interfering proper- 
ties. We will assume that the materials are sufficiently similar in type to avoid 
difficulties of this sort. But we do allow for a functional relationship between the 
within-cell standard deviation for a given material and the magnitude of the 
measurement for this material. 

Such a relationship is best detected by ordering the columns according to 
increasing order of magnitude of their averages and studying the within-cell 
standard deviation for trends from left to right. If such a trend is found to exist, 
the pooled within-cell standard deviation for each column may be plotted 
against the average of all measurements in this column, and a simple curve 
fitted to the points thus obtained. Generally, a straight line will be sufficient. 
Thus, denoting by o,, the standard deviation within cells in column j and by 
9; the average of all observations in column j, one will have the approximate 
relation: 


— I 
DL Yin — Gis)” 


= = : : in i=] k=1 
(2) o., = A + Bg; where o,, is estimated by ¢.,; ‘man-D 


By making a transformation of each observation Yiiz into 2;;, such that 
{3) z = K log, (A + By) —G 


where K and G are arbitrary constants, chosen for computational convenience, 
the trend expressed by equation (2) will be largely eliminated. Indeed, if « is 
the within-cell error for the transformed variable z, it follows from (3) by the 
common rules for the propagation of error that a good approximation for o, 
is given by: 


.= _— = ex anienieaeeeleks. gies KB 


It will be assumed in all that follows that a transformation of the type given 
by equation (3), or any other convenient type, has been applied whenever 
necessary, so that for the transformed variable, the standard deviation within- 
cells may be considered to be essentially independent of the magnitude of the 
(transformed) measurement. According to our assumptions, it is then constant 
over all cells of Table I. The mathematical model and the analysis of the data 
is to be thought of in terms of the transformed variable z;;, . 
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THE MATHEMATICAL MopEL 


Our aim is to express 2;;, in accordance with a relation of the type given by 
equation (1), where M is now the transformed measurement z;,, . Keeping in 
mind that the values of Q are assumed not to be known, any useful model must 
be based on certain assumptions restricting the generality of equation (1) and 
allowing the variable Q to be eliminated. It must also be borne in mind that the 
variance of those variables in C that are responsible for within laboratory 
variability has been stabilized with respect to the magnitude of the measure- 
ment by means of a suitable transformation. 

The basic assumption made in this paper is that, apart from within laboratory 
error, the results obtained by the various “laboratories” are linearly related to 
each other. 

We will first consider the case in which no I-factors are present, and extend 
the model subsequently to the general case. We have represented by ¢;;, the 
within-cell error corresponding to the transformed measurement z,;, , and 
assume that the population mean of ¢,;, is zero. Denote by ¢;, the expected 
value, with respect to ¢;;, , of the observations in the (2, 7)th cell. Let £; be the 
average of the [;; in column j and represent by yu, the average of the ¢;; in row 7. 
Thus: 


(5) i= ¥ tase 


i=1 


(6) i= 2 $:s/b 


Furthermore, let c be the average of all ¢;; in Table I: 


a 


(7) c= 2 Di tu/ab=F=a 


t=1 j 


where ~ and 7 are the average values of the £; and the u;, respectively. 

From the assumption of linear relationships among the a laboratories it fol- 
lows that the values obtained by each laboratory are linearly related to the 
corresponding average values of all laboratories. Thus, there exist a constants 
vy; and a constants 6;, such that: 


(8) fi = 15 + BE; — 


Summing (8) over all j yields: 


De fii = bys = by, + 8; » & -—) = by, 


hence: (9) vy; = w;. Summing (8) over all 7, and using (7) and (9), one obtains: 


Do fii = ag; = ac + & — 0) » 8: 


Now, since the materials are not all alike with respect to the measured property, 
the quantity &; — c (or &; — 2) cannot vanish for all"j; hence the preceding 
equation yields: 
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(10) 2 B; = a, B 


By definition of ¢;; we have: 

(11a) Zin = $5 1 Cczn 
Introducing (8) and (9) into (11a) yields: 

(12a) Zine = wi + BE: —O + Gir 


Equation (12a) represents the model in the absence of I-factor effects. If 
I-factor effects are present, a modification of the model becomes necessary. 
Among the interfering properties there will be some that affect the measuring 
process in identically the same manner as the property of interest, Q. Since our 
model does not involve absolute standards for Q, interfering properties of this 
type will be indistinguishable from Q and go unnoticed in the study. On the 
other hand, some of the interfering properties may affect the measurement in a 
manner different from Q, so that different laboratories will respond to their 
presence in a manner unlike their response to Q. 

Let \;; represent the effect of the totality of interfering properties of the 
latter type for material j, on laboratory 7. In order to make allowance for the 
effect of \ in our model it will be assumed that for any column j the material in 
that column is a random selection from a large population of materials all 
having the same value of Q, but representing many values of the interfering 
properties, such that the average of \,;; over this population is zero in each 7, 
j cell.’ 

The generalized model is obtained by substituting the following equation for 
equation (lla): 


(11b) Ze = Sig HYG H+ Ge 

which leads to: 

(12b) Zizzx = wi + BE: — OC + AG + Eze 
Representing the cell-average of z,;, by z;; we can write: 


(13) 2:3 = wi + BME; —0C) +A; + dX €:54/N 


The random variables \,;; and >>, €;;,/n may be grouped together. Let: 
(14) ni = AG + ke €;;./n 
k 


If the symbol V(u) represents the variance of u, we have: 


(15) Vinii) = Vii) + Vein)/n 


1 Strictly speaking, the random effect \;; may also include a component unrelated to the 
materials and associated with the equipment or the procedure. An example of this is found 
in calibration errors of an instrument that are randomly distributed along the range of readings 
of the instrument. No misinterpretation of the analysis will result from the presence of such 
components, provided that this possibility is kept in nnd. 
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The transformation of y into z has resulted in the stabilization of the variance 
of ¢;;, . We will assume furthermore, that V(A;;) is likewise independent of 7 
and j. Thus, (15) becomes: 


(16) Vin) = VO) + Vn 
From (13) and (14) we obtain: 
(17) 2:5 = wi + BME; — OC) + 05; 


ANALYSIS OF THE DaTA 





One of the main purposes of the study of a measurement process is to evaluate 
its precision under practical circumstances. The cost and effort of an inter- 
laboratory study are justified only if it leads to an evaluation of the agreement 
that may be expected between results obtained in different testing institutions. 
But, in view of our interpretation of experimental error, such an evaluation is 
not only of practical value, but also leads to a better understanding of the 
capabilities and limitations of the measuring process. 

In the preceding sections of this paper we have established that at least 
three sources are responsible for possible discrepancies between results: 


1. The error expressed by V(e), characterizing the variability among replicate 
measurements. 

2. The error expressed by V(A), characterizing the differential response of 
different laboratories to interfering properties. 

3. The error due to systematic differences between laboratories characterized 


in our model by the variation of two parameters, yu; and 8; , i.e., by V(u) and 
V(8). 


The main problem is to evaluate the relative contributions of these various 
sources of error to the overall variability of the results obtained for any given 
material in different laboratories, and to estimate this overall variability. It 
should be borne in mind that the breakdown of the overall variance is not 
necessarily the same over the entire range of values taken by the various ma- 
terials. It is for this reason that the analysis must be carried out first for a 
fixed value of j (a given material) and then studied in terms of varying j-values. 

We will use the following notation and terminology. Given any finite set of 


numbers u,; , U2, °** , U, , We define their variance as: 
(18) Vw) = > (u; — a)?/(n — 1) where @ = * u;/n 
We define the co-variance of a finite set of paired quantities (u; , v1); (ue , v2); 
s+ 5 (uy, t,) a8 
(19) Cov (u,v) = 7 (u; — av; — d)/(n — 1) 


Given any function of the z;; , 


f(z;3) = fle: a's Ey »Zab) 
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we consider it as a function of the fixed quantities ¢;; and the random variables, 
of zero expectation, 7:; : 


fess) = PGi, +++ 2a 8+» Zan) 


= f(Su H+ mr, ce? Ses HF es ot? Sea H Mes) 


By the “expected value” of f(z,;;), denoted E[f(z,;)], we mean its expected 
value with respect to the random variables 


ms * OS gE Sg ee 
It is easily shown that if f(z;;) is quadratic in all z,; , 
(20) Effi] = fii) + Elf(nis)] 


Now, consider equation (12b) and let 7 have a fixed value. Then it is easily 
proved that: 


(21) E[V(@.sx)] = Vu) + & — °V(B) + 2G; — ¢) Cov (u, 6) + VA) + VEO) 


Estimates for V(x), V(8), &; and ¢ are obtained as follows: 
Let the row and column averages of the z;; (Table I) be represented by m; and 
x; , respectively: 


(22) mm; = Dib t= D 2;/a 


Evidently, m; is an estimate of »,; , and x; is an estimate of ¢; . Furthermore, 
the quantity m = = )>; >); z;,/ab is, of course, an estimate for c. We now 
treat the z;; data in accordance with the usual analysis of variance technique 
for two-way classifications, but interpret the mean squares in accordance with 
the model expressed by equation (17) together with conditions (7) and (10). 

Using equation (20) and equations (7), (10), and (17) the expected values in 
Table II are readily derived. 

Frem these equations it is possible to derive estimates for V(u), V(@) and 
V (é), provided that an estimate of V(n) is available. Now, the analysis of variance 
does not provide an estimate of V(n), but such an estimate can be obtained by 
applying the theory of linear regression to equation (17). This will be done 
below, and it will also provide an estimate for V(A), by using equation (16), 
since V(e) can be estimated by means of equation (4). 


Taste II 


Degrees of 
Source of Variation Freedom E (Mean-Square) 


Laboratories a-1l bV(u) + V(n) 
Materials b-1 aV(é) + V(n) 
Interaction* (a — 1b — 1) V(B)V(E) + Vin) 


* Laboratories X Materials 
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The covariance term in equation (21) makes allowance for a possible relation 
between the two parameters yu; and B, of the response line for laboratory i. 
From a physical point of view such a relation may well be expected in many 
cases. Suppose indeed that the measuring process includes a “single-point 
calibration”, i.e., a calibration in terms of some standard weight or sample at 
one given point in the range of Q values. Then all laboratories must agree at 
this point to within V(7). This in turn implies that the response-lines for all 
laboratories pass through a common point, i.e., that they are concurrent. It is 


readily seen that such concurrence implies a linear relationship between the 
parameters py; and 8; . 


B: = alu -O +k 
Since >>; 8; = a, it follows that K = 1, hence: 
(23) 8B; =a; -)O +1 
Introducing this expression for 8; into the model (17), we obtain: 
25 =a tatu -—O8; —-O +E -—O + 1; 
which can be written: 


(24) 2i=ecttu -O+(; —O tole; — OF; —O + 13 


An analysis of the model expressed by equation (24) has been given by Ward 
and Dick [2, 6] who propose an iterative method for the estimation of a and of 
the variance of the error term 7,;; . The-first step in the iteration is identical 
with a procedure first proposed by Tukey [5] which consists in extracting from 
the interaction sum of squares in Table II a portion corresponding to one degree 
of freedom. This single degree of freedom measures the significance of the 
multiplicative term a(u, — c) (€; — c). In practice the first step in the iteration 
is sufficient. The estimate of the parameter a by Tukey’s method is given by 


d De 2i(m, — m)(x; — 2) 


a= Lm, — ty 2, — 2 


and the corresponding sum of squares by 


ia De 2ii(m; — m)(x; — #)) 

a (m; — my)’ a (x, — 4)" 
Frequently the data indicate a tendency of all laboratories to show better 
agreement at one point of the range of values without, however, coming close 


enough to justify the model of strict concurrency of all laboratory-response 
lines. Such a situation can be represented by the equation: 


(25) B; = alu — c) + K; 


which differs from (23) by allowing K, to vary from laboratory to laboratory. 
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From i 8; = a it follows that Du K, = a. Let 


(26) ay. 1 + 6; 
Then >>; 5; = 0 and the model expressed by equation (17) becomes: 


(27) 2;=ct+tw—-9O+€6 -—O tatu; —O€ — oO + [n; + 8; — 0] 


Applying Tukey’s method to this model, one will obtain an error estimate 
whose expected value will exceed that of V(), because of the variability due 
to the 6; . Consequently, a comparison of this error-estimate with the estimate 
for V(n), to be obtained below, will be a clue to the presence of 6; terms differing 
from laboratory to laboratory. 

Using equation (27), the breakdown of the variance of z;;, , for a given j, can 
be expressed in a more useful form than that of equation (21): 


(28) = E[V@in)] = [1 + a&; — OPV) + & — 0°V(8) + VO) + VEO 


In this equation, the effect of V(@) has been replaced by two quantities, a and 
V(6), providing a more sensitive method of analyzing the total variance. 

An estimate for V(é) could be obtained by comparing the error-estimate 
provided by Tukey’s procedure with the estimate of V(n). More conveniently, 
V(6) can be estimated by means of the following relation, which follows from 
(25) and (26): 


(29) V(8) = a’ V(u) + V(6) 


This equation is solved for V(6), using the previously obtained estimates for 
V(8) and V(u), and the estimate of a provided by Tukey’s procedure. 

There remains the task of finding an estimate for V(n). 

From (22) it follows that: 


(30) 


Consequently: 


ee Nii eo zs Nii 
(31) €; —oc =(;-% -|- ae a 


a ab 


Equations (17) and (31) establish, for any fixed value of 7, a functional linear 
relationship between z;; and x; . Both these variables are subject to error and 
their errors are correlated. Introduce the variable: 


(32) ti; = 02; — 25; 


Combining equations (17), (31), and (32), it is seen that for any fixed value of 
i, the variables z;; and #;; are linearly related, and that their co-variance with 
regard to the errors 7;; is zero. The variance of the error of z,;; is V(), while 
that of the error of ¢;; is (a — 1)V(»). Using the method described by Lindley 
[3] for the estimation of a linear relationship between two variables subject to 
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error, with zero correlation between their errors, one obtains, for any given 7, 
the following estimates. 


Yi 
(33) f+ Oe, 


(34) P(n) = = +5 ld + v2 — v.AQP)] 


where 


2 1 
(35) 4 “8 Fir 


_ oX — 2aP; — (a — 2/4, 
(36) = "Fa — NeP, — 2) 


(37) X= VG; - 9’ 
(38) Z; = so (i; — m,)” 
(39) P; = a (2,; — m,)(x; — 2) 


Equation (34) provides an estimate of V(n) for each value of 7. These values 
are averaged to provide the overall estimate of V(n). The V(n) estimate obtained 
by this procedure is used for two purposes. It is introduced in the expressions 
given in Table II, which now suffice to obtain estimates for V(u), V(8), and 
V(é) (this last parameter is of no relevance to our problem). The estimate of 
V(n) also serves to derive an estimate for V(A), using equation (16) and the 
estimate of V(e) given by an equation of the type of (4). 

We are now in a position to estimate the overall variance V(z,;.) for each 
material j, by substituting in equation (28) the estimates for all the quantities 
occurring in the second member. (The quantity (, — c) is estimated by (x; — Z)). 

The relative contribution of each of the component parts to the total variance 
varies in general with the magnitude of the measured value for the material 
under study. It is useful to tabulate these relative components of the total 
variance for various ranges on the scale of measurement. If it is desired to 
express the variability in terms of the original scale, y, this is easily accomplished 
by remembering that 


ieee 
A + By 
1 
72 


K E + | V2) 


By means of this relation, the total variance is expressed in the original 
scale. The ratio of the variance due to each component to the total variance 
is unaffected by the transformation of scale. 


og, = K 


Cy 


hence: 


(40) V(y) = 
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THE PERFORMANCE OF INDIVIDUAL LABORATORIES 


It is often of interest to compare the results obtained by a given laboratory 
with the average result for all laboratories. Equation (33) provides the slope for 
the ith laboratory. The quantity m; , equation (22), is an estimate of yu; . Finally, 
equation (34) is an estimate of V(n), based essentially on the data of the 7th 
laboratory. Thus, the two parameters defining the response line for the ith 
laboratory, as well as the scatter about this line, are available. 

It should be pointed out that the values u; and 8; may not be stable with 
respect to time. A repeat experiment, carried -out at a later date, may well yield 
appreciably different values for u; and 8; for any given laboratory. Only experi- 
mentation can provide information of this sort. It is, therefore, important to 
follow up on the results provided by the initial interlaboratory study, by means 
of subsequent studies on a smaller scale, using perhaps two or three materials 
only. 


THE IMPROVEMENT OF PRECISION 


If a series of tests has shown that the values of yw; and 6; are stable for each 
laboratory, then the results obtained by any given laboratory for any material, 
j, can easily be “adjusted” to conform more closely to the average of those 
obtained by all others by solving the equation 
(41) Zin = A + BE; — é) 
for &; , using the observed value of z;;, - 

Except for the random fluctuations in the estimates f; and 8; , the results 
obtained by the various laboratories will then agree to within V(A) + V(e). 
Further improvement is possible by averaging replicate measurements, In this 
manner V(e) is divided by the number of replications, and it is seen that the 
ultimate precision attainable is given by V(A). 

Should the values of yu; and 6; , for any given laboratory, fail to remain con- 
stant with time, then they cannot, of course, be used in the manner just outlined. 
However, assuming that during any given testing session in a particular labora- 
tory there prevails a given line, with parameters u; and 6; , these parameters 
can be estimated during the same session by means of two “reference samples”, 
i.e., materials for which the values £; are either known or have been assigned 
definite values. The practice of correcting experimental results by means of a 
single reference material is quite common. The theory developed in this paper 
shows that in the general case, it is necessary to resort to two reference materials 
in order to obtain the best possible adjustment of the data. However, if the 
analysis has shown that the lines corresponding to all the laboratories are con- 
current, it will suffice to resort to a single reference material which should be 


chosen in such a way that its value is as far as possible from the point of con- 
currence. 


A NuMERICAL EXAMPLE 


As mentioned earlier, the practical aspects of the theory developed in this 
paper are dealt with in a companion paper [4]. Nevertheless, in order to give 
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Taste III 
SULFUR IN PETROLEUM* 


LABORATORY MATERIAL 


0264 
.0282 
-0271 


.0286 . .189 
.0274 . P 
.0276 .. . 187 
Average(1) .0107 .0282 . .189 . .912 
Average Range(1) .00075 .00270 . -00642 .0165 .0235 


(1) Of triplicates within cells. 
* Reference: C. C. Hale, E. R. Quiram, J. E. McDaniel, and R. F. Stringer, 
“Eeso Lamp Method for Sulfur’, Anal. Chem. 29, p. 386, (1957). 


the reader a more concrete picture of the proposed method, an example is 
appended. 

Table III is taken from a study of a method for the determination of sulfur 
in petroleum [Anal. Chem., 29 p. 386 (1957)]. Triplicate determinations were 
made for six materials by each of four laboratories. The empirical relation 
between the average of triplicates and their standard deviation* inferred from 
these data is ¢ = 0.000252 + 0.02565 y. The transformation used was: 


z = 10°[4 + log, (0.000252 + 0.02565y)] 


The analysis of variance of the transformed data is given in Table IV, and the 
results of the regression analysis in Table V. The following estimates were 
derived from the analysis: 


Vi) = 124.4 V(8) 
Vr) = 66.4 Viu) 
& = —1.070 x 10° 


* In the example here considered the estimates of the standard deviations were obtained 
from the average ranges, rather than by the usual sum of squares formula. 
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TABLE IV 
ANALYSIS OF VARIANCE OF SULFUR DATA” 


Source of Degrees of Mean 
Variation Freedom Square 


Materials 5 1,470,018 

Laboratories 3 209. 

Materials & Laboratories 15 12.1 
Term in a 1 220.3 
Remainder 14 04.3 


“) In transformed scale 
2) See eq. (24) and references (2, 5) 


Finally, in Table VI, the variance of z is broken down into four components 
in accordance with equation (28), where among the conclusions of practical 
usefulness that can be drawn from this analysis, the following two are probably 
the most important. 

1. The ratio of V(e) to V(A) is approximately equal to two. Consequently, 
three or more replicates suffice to render the contribution of the e-error less than 
that of A. Since V(A) cannot be reduced by replication, a number of replicates 
considerably larger than three is therefore of no practical value for the improve- 
ment of precision. : 

2. The contribution of the between-laboratory variability to the total vari- 
ability (based on single determinations per material per laboratory) is negligible 
for materials containing more than 0.1% sulfur. For materials containing 
between 0.01 and 0.1% sulfur, the between-laboratory variability is more pro- 
nounced. If it were desirable to improve the- overall precision of this method, 
especially in the range of low sulfur contents, this could be accomplished by 
means of a reference sample whose sulfur content is of the order of 0.01 per cent. 
Specifically, each laboratory would calibrate its results versus the reference 
sample by plotting, on the transformed scale, its result for the reference sample 
versus the assigned value, and drawing a straight line through this point and 


TABLE V 
LINEAR REGRESSION ANALYSIS OF SULFUR DATA” 


Laboratory a B v% n) 


1549.3 .9893 48.7 
1535.3 1.0009 237.3 
1539.1 1.0097 118.9 
1541.1 .9999 26.4 


Average 1541.2 1.0000 107.8 


() In transformed scale. 
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TABLE VI 
RELATIVE COMPONENTS OF VARIABILITY FOR SULFUR DATA 


Percent of total variability due to 


% Coefficient 
Material V(4) V(u) % Sulfur of Variation“ 


KSESSS 


34.8 


“® See equation (28) 
(2) Based on a single determination per laboratory and including all sources of 
variation. 


a point on the bisector of the axes corresponding to 1.16% sulfur (the point in 
which the lines of all laboratories already concur. It is readily shown that for 
V(é) = 0, the point of concurrence has an abscissa and an ordinate both equal 
to c — 1/a). The corrected value corresponding to any value found by this 
laboratory is the abscissa of that point on the straight line for which the ordinate 
is the observed value. 
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Factorial Experiments in Life Testing’ 


Marvin ZELEN 
National Bureau of Standards 
Washington, D. C. 


This paper discusses procedures for analyzing factorial experiments, where the 
experiment deals with the life testing of components or equipment. These procedures 
assume an underlying general distribution of ‘“‘times-to-failure’’, of which the ex- 
ponential, Weibull, and extreme value distributions are special cases. Statistical tests 
and confidence procedures are outlined, and an example illustrating the procedure 
for life-test results of glass capacitors is included. Small sample approximations, 
which are adequate for practical applications, are given for the proposed procedures. 
This is shown empirically by generating thousands of life-test experiments on an 
electronic computer. An empirical sampling investigation is given of the robustness 
of the proposed procedures. From the sampling results, it is concluded that these 
techniques are sensitive (non-robust) to departures from the original assumptions 
on the probability distribution of failure-times. An investigation is also given of a 
transformation which appears to give robust results. These same techniques carry 
over exactly to the situation where one is analyzing an array of variance estimates 
from an underlying normal population. 


1. INTRODUCTION 


Recently much attention has been devoted to the study of life-tests. These 
so-called life-test experiments are characterized by running one or more nominally 
identical components (or pieces of equipment or parts) under a given set of 
operating conditions and noting the number of hours of satisfactory performance 
for each component. In the case of electronic equipment the data recorded will 
usually be the “‘time to failure”, where “failure” is defined as that state of the 
item where satisfactory performance is no longer achieved. In the case of fatigue 
testing, such as tests on ball bearings, the data recorded will be the number of 
revolutions until failure, where a failure is a malfunctioning of the bearing. 

Many of the methods and techniques associated with life-tests assume that 
the distribution of failures follow the exponential distribution, cf. Epstein and 
Sobel [8], [9]. These techniques, however, apply only to situations where one 
is testing all components under the same fixed environmental conditions. The 
object of this paper is to present methods for analyzing failure data when the 
data are taken over a range of different environmental conditions. In engineer- 
ing parlance these tests are often referred to as ‘‘combined environment” tests. 

The methods presented here deal with the analysis of life-test data where the 
data conform to a factorial experiment. These techniques will enable one to 


1 This paper is based in part on work done while the author was at the University of Cali- 
fornia (Berkeley), Spring 1957. 
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evaluate the effect of different environmental conditions on the failures. For 
example, it is a common rule of thumb to state that the “life” of components 
made chiefly from inorganic materials (such as some capacitors) is halved for 
every 10°C rise in temperature. This paper will show how such rules of thumb 
can be checked against empirical data. 

A not uncommon procedure in life-testing is to accelerate the failure of com- 
ponents by testing at an extreme environmental condition.(such as high temper- 
ature) and then attempt to extrapolate the accelerated test data to regular 
operating conditions. This paper will present criteria which must be satisfied 
if such extrapolations are to be meaningful. 

The techniques and methods presented here were originally motivated by 
life-test problems. However these same formulas apply to the analysis of a 
factorial array of variance estimates when the underlying observations follow 
normal distributions. 


2. OUTLINE OF EXPERIMENT 


Often the experimenter desires to investigate the effects of several factors 
on the performance of a piece of equipment or on the behavior of components. 
It is widely accepted in many experimental sciences that a considerable ad- 
vantage is gained if such multi-factor experiments are conducted so that the 
effects of changing each factor can be evaluated jointly with the effects of the 
other variables which might conceivably affect the outcome of the experiment. 
Life-test experiments are no exception. 

For the purpose of exposition, this paper will only discuss experiments in- 
volving two factors, however the extension to more than two factors is straight- 
forward. Let the two factors be denoted by A and B and the number of levels 
(or different conditions) for each factor be a and b respectively. Then a treatment 
combination can be denoted by the symol (77) where 7 and j refer to particular 
levels of factors A and B respectively. There are a levels of factor A, i = 1, 
2, --- , a; similarly for factor B we have j = 1, 2, --- , b. Thus the total number 
of different treatment combinations is ab. For example, factor A may refer to 
voltage and factor B to temperature. If there were 4 voltages under investi- 
gation, and 2 different temperatures, then a = 4, b = 2, and there would be 
ab = 8 different treatment combinations. 

The experiment will be planned so that for each treatment combination, n 
components are simultaneously placed on test. The tests for each treatment 
combination will be terminated when exactly r < n of the test items have failed, 
where r is a number chosen in advance of the experiment. Since the first failure 
will correspond to the smallest failure time, the second to the next smallest 
failure time, etc., we will denote the r failure times for each treatment combi- 
nation byt; << +: <t.(r <n). 


3. ASSUMPTIONS UNDERLYING FAILURE DISTRIBUTIONS AND 
OBJECTIVES OF THE ANALYSIS 


In many studies of failure data, the exponential distribution has been used 
to describe the distribution of failures. We shall initially assume that the failure 
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times follow the exponential distribution having the probability density function 


1 -«-eye 
1) f(t) = 9° for t>G 


0 for t<@: 


The parameter @ is usually called the ‘“‘mean-time-between-failures”’ ; its reciprocal 
6~* is termed the “failure rate”. The parameter G is referred to as a “guarantee 
time” and has the property that ne failures occur before time G is reached. In 
many applications G is set equal to zero (i.e., G = 0) resulting in a one parameter 
exponential distribution. We shall assume that the underlying distribution of 
failures for the (7j) treatment combination is exponential with parameters 6;; 
and G;; . The object of the analysis will be to analyze the effects of factors A 
and B on the mean-time-between-failures, (0;;), irrespective of the value of the 
guarantee time (G;;). 
For this purpose we shall write 0;; as 


2) 1; = ma,b¢;; 


where the quantities a; , b; , and c;; are positive constants and subject to the 
restrictions 


a b a b 


3) Ila = [1] = Ile; = [Te, =1. 


i=1 i=1 i=1 i=1 


These parameters are similar to the main effect and interaction terms one as- 
sociates with the usual analysis of variance models. Note that log 6;; = log m + 
log a; + log b; + log c,; which is exactly the analysis of variance model for 
a two-way classification. The quantity m is a constant common to all treat- 
ment combinations; a; represents the contribution from level 7 of factor A; b; 
is the contribution from level j of factor B; and ¢;; is a constant which depends 
on both factors. Thus we can regard the a; and b; as main effect terms and the 
c;; aS an interaction. The restrictions in equation (3) are imposed so that the ab 
different @;;’s can be expressed as a function of exactly ab other independent 
parameters. 

Our procedure for analyzing the failure data will consist of (7) estimating the 
values of the parameters m, a; b; and c;; and (ii) making appropriate statistical 
tests of significance or confidence statements on the parameters of the failure 


distribution. Some of the pertinent hypotheses and their implications are given 
below: 


Hypothesis Model Implication 


a,b;c;; = 1 6; =m 6 not dependent on different levels of factors A and B. 
bje,; = 1 6 only depends on the levels of factor A. 
a,¢;;= 1 6;; = mb; 6 only depends on the levels of factor B. 
¢;=1 6; The effects of factor A can be evaluated inde- 
pendently of the effects of factor B and vice-versa. 
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When 6,; = ma,b; , we will say that the effect of the two factors is completely 
multiplicative. With this model, accelerated environment testing seems feasible. 
For example, factor B might refer to temperature with one of its levels chosen 
to be an unusually high temperature in order to increase the failure rate. Then 
we will still have information on the parameter a; which will apply to operations 
under regular conditions. Thus a desirable condition for accelerated environ- 
mental testing is that the model be completely multiplicative. 


4, FoRMULAS FoR EstIMATEs, TESTS OF SIGNIFICANCE, 
AND CONFIDENCE STATEMENTS 


4.1 Estimation formulas 
With respect to any treatment combination there will be r (r < n) ordered 
failure times ?#,; < tp < --- < #,. Define 


4) g — Mi— 4) + M—  — 4) 


r—l 


i = > Re. 


For the (ij) treatment combination, let the 6 quantity be denoted by 6,;; . Then 
the (maximum likelihood) estimates for the various parameters are given by 
(5). The quantity 6;; is actually the minimum variance unbiased estimate of 0;; . 


mh = (TI II a) 


t=1 j=1 


A caret (~) is used to indicate that the above quantities are estimated from 
the data. Note here we compute geometric averages rather than the usual arith- 


metic averages. In order to compute these quantities, it is best to use logarithms, 
i.e., 


a b 

log * = 5 > >d log 6,; , ete. 
ab i j=1 

Hence the log quantities are similar to the estimates of the grand mean, main 

effects and interactions for the usual two-factor factorial experiment. The 

extension to more than two-factors proceeds in this same manner. 
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4.2 Tests of significance and confidence statements 


The formulas for making appropriate tests of significance are summarized 
in Table 1. The test criteria for the various hypotheses all take the form 
M = A/C. We will term M the modified likelihood ratio. The quantity M is 
approximately distributed as a x’ variable with the indicated degrees of freedom. 
This approximation seems suitable for practical use. Large values of M indicate 
significant departures from the null hypothesis. Better approximations are 
discussed in section 6.4. 


TaBLe 1. Formulas for Making Tests of Significance!.?.4 


Degrees of 
Null Hypothesis freedom Ci 


3ab — 2a +1 
Ay = 2ab(r — 1) log . ere a 
3ab — 2b +1 


Az = 2ab(r — 1) oe ( Crm1l+ 6ab(r — 1) 


ab + 4a + 4b — 11 
oe ai , re = aa , oa enna nainaoatdeagalines 
(a —1)(6 —1) As = 2ab(r — 1) oe ( 2 Cs =1 + Gabe — 1) 


1 1 tt. 
A, = 2ab(r — 1) - Dy log (3 Zz bas) Ci = 1+ aon - 
‘ j - 


1 1 a+l1 
:aiciy = As = 2ab(r —1)- 1 - i035 = ee 
aicij = 1 5 = 2ab(r 14 St (2 ¥ aes) Cs =1+ 3 


1 ab +1 
: asbjcig = 1 Ag = 2ab(r — 1) 1 - : 
asbjcig 6 ab(r ) log (3 E abs) Ce =1+ <5 


1. All logarithms are to base e. 
2. A; is likelihood ratio test and is asymptotically distributed as x*. 
3. Mi = A,/C; is modified likelihood ratio test and is approximately distributed as x?. 


With respect to applications, the hypotheses that 0;; = ma,c,;; or 0;; = mb,¢;; 
seems of doubtful practical value. However for completeness these formulas 
are also included in Table 1. 

Frequently one desires to compute confidence intervals for the ratio of the 
main effects a;/ai (¢ + 2’) or b;/b; (j + 7’). Tables 2 and 3 can be used for this 
purpose. Tables 2 and 3 give the appropriate values for 99% and 95% confidence 
intervals respectively. Each of these tables depends on (i) the number of 6 
terms (g) in the geometric average of 4, or 6; , and (ii) the quantity » = 2 (r — 1). 
The limits for the two -sided confidence interval for a,;/ai take the form 

me 1 (4; 
lower limit: Z ( 4) 
4, 


upper limit: Z (4) 
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where Z is obtained from Table 2 or 3. The number of terms in the geometric 
average 4; isg = b. When placing confidence intervals on b;/b/ , we have g = a. 


TABLE 2 
Constants Z for Two Sided 99% Confidence Limits on a;/a/, or b;/b} 


Number of terms Degrees of freedom vy = 2(r — 1) 
in geometric 
average (9) 


5.89 
4.59 
3.89 
3.44 


5.151 
For larger values of g and »v use Z = ex ee 


p 
Vo» — 1) 


TABLE 3 
Constants Z for Two Sided 95% Confidence Limits on a;/a{ or b;/b/ 


Number of terms Degrees of freedom »v = 2(r — 1) 
in geometric 
average (g) 


3.48 
2.77 
2.40 
2.19 
2.05 


3.9200 
V ad» — 1) 


For larger values of g and »y use Z = exp 


5. NUMERICAL EXAMPLE 


Life-test data were obtained for glass capacitors over a range of different 
temperature and applied voltage conditions. The environmental conditions of 
the experiment were: 


Factor A: applied voltage: 200, 250, 300, 350 volts 
Factor B: temperature: 170°C, 180°C. 


Eight capacitors were tested for each of the eight voltage-temperature combi- 
nations. After the fourth failure, the tests were suspended for every treatment 
combination. Thus we have a = 4, b = 2,n = 8, r = 4. Table 4A summarizes 
the data. 
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TABLE 4A 


Life Test Results of Capacitors (Hours) 


Temperature 


170°C 


Computational steps for analysis 


200 


Applied Voltage 


250 


572 
690 
904 
1090 


315 
455 
473 


300 


315 
439 


i) For each voltage-temperature combination compute 6,;; using equation (4). 
Table 4B summarizes the results of these calculations. 

ii) Compute geometric averages, using the 6’s, for each temperature and 
voltage. Also compute the grand geometric average. These calculations are also 


included in Table 4B. 


ii) Using the formulas of equation (5) with n = 8 and r = 4, compute m, 4; , 
6; , and é,; . The result of these computations are: 


200 
50 


6,, = (mh = 556] x 7 
300 


350 


200 
x 250 
300 


350 


Voltage 
d, = 1 .1680 


ad, = 1.3329 
ds; = 0.7229 
ay = 0 .8890 


170°C 
1.5173 


Ey; = 0.9087 
és, = 0.9309 
és, = 0.7793 


a —_— 
“= 


Temperature 


‘ 170°C ‘s = 1.5057 


180°C | 6, = 0.6641 


180°C 
= 0.6589 


1.0999 
1.0746 
= 1.2843 


Thus if we desired to predict the value of 6;; for 250 volts and 170°C temperature, 


we would have 


dzb,é., = (556)(1.3329)(1.5057)(0.9087) = 1012.2 hours. 


This answer should agree with 6,, except for rounding errors. 
iv) Finally using the formulas in Table 1, we compute the various tests of 


significance. These are summarized in Table 4C. 





Temperature 
170°C 
180°C 


Geometric 
Average 


Hypothesis 


6:3 = max; 
6:3 = ma; 
6:3 = mb; 
6: =m 
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TasLe 4B 
Summary of 6;; 


Applied Voltage 
Geometric 
250 300 Average 


1013 563 837 
541 287 369 


741 402 556 


TABLE 4C 
Summary of Tests of Significance 


x? (critical 
value at .05 level) 


7.82 
9.45 
12.59 
14.07 


Note that the hypothesis 6;; = ma,b; is supported by the data (i.e., M = 
1.23 < 7.82) which indicates that the model may be completely multiplicative 
over the range of conditions of this experiment. 

v) For this example it is of interest to calculate the change in the average 
failure rate for a 10°C rise in temperature. The ratio of the average failure 


rates is 


and indicate that the decrease of 10°C resulted in a failure rate .441 times that 
for 180°C. Since this ratio is a function of the life data and is subject to variability, 
it would be desirable to have a confidence interval for this ratio. Adopting a 
95 per cent confidence coefficient, we can enter Table 3 with g = a = 4 and 
v = 2(r — 1) = 6, and find Z = 2.41. Hence the 95 per cent confidence limits 
for the ratio of the failure rates (b./b,) can be obtained using equation (6). 


This results in 


_1 


lower limit: 24l 


(.441) = .18 


upper limit: (2.41)(.441) = 1.06 
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6. MISCELLANEOUS REMARKS 


6.1 Guarantee time equal to zero 


In some applications the guarantee time G is taken equal to zero. This requires 
a slight modification of 6 (equation 4), the tests of significance (Table 1), and 
the tables for placing confidence intervals on ratios of main effect terms. Table 
5 summarizes these changes. 


TABLE 5 
Summary of Changes in Formulas when G = 0 


Equation or Table Number Change 


Equation 4 Define § = [7 + (n — r)t,]/r 
Table 1 Change (r — 1) tor for the A; 
and C; columns 
Tables 2 and 3 Change vy = 2(r — 1) tov = 2r 


6.2. Application of formulas to analysis of variance estimates 


Occasionally the situation may arise where one has a factorial experiment 
such that there exists an independent variance estimate for each treatment 
combination. It is desired to analyze the effects of the factors on the variances. 
The formulas of the preceding sections can be used for this purpose provided 
the variance estimate s7; for the (ij) treatment combination follows a 07; x’/v 
distribution with » degrees of freedom. Then the model analogous to equation 
(2) is 


8) O45 = ma,b,c;; . 


Some of the formulas have to be changed slightly. These changes are summarized 
in Table 6." 


TABLE 6 
Summary of Formula Modifications when Analyzing Variance Estimates 


Equation or Table Number Change 


Equation 5 Use s?; instead of §;; . 
Table 1 i) Column A; : change 2ab(r — 1) 
to aby 
ii) Column C; : change 6(r — 1) to 
3v. 
Tables 2 and 3 Use the degrees of freedom v 
for the column entry. 


1In connection with analyzing sample variances, Professor Robert Bechhofer has also 
considered this same multiplicative model. 
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6.3 Applications to other than exponential distributions 


The techniques described in the first five sections of this paper were predicated 
on the assumption that the failure times followed the exponential distribution. 
The question naturally arises as to what procedures should be used when the 
underlying distribution of failure times is other than exponential. For example, 
the Weibull distribution has been widely used in connection with describing 
the failure distribution of ball bearings, cf. Lundberg and Palmgren [13] and 
Lieblein and Zelen [12]. Kao [11] has also used the Weibull distribution to 
empirically fit the life distribution of certain vacuum tubes. 

Fortunately, all results for the exponential distribution can be used for 
distributions having the probability density function of the form 


9) HH JO oxy -[ h(2) ac} for t>0 
0 for ¢ <0. 


The function h(é) is simply a non-negative function of the time to failure ¢. The 
exponential and Weibull distributions are special cases of the general family of 
distributions described by equation (9). 

The quantity h(¢) is sometimes referred to as the “hazard”, “instantaneous 
rate of failure’’, or ‘conditional failure rate’’. Its interpretation is that if a com- 
ponent is put on test at time zero and performs satisfactorily up to time ¢ (say), 
one might ask about the probability of the component failing in the next instant 
(say) between time ¢ and ¢ + At. This probability is the conditional probability 
of failure and implies that the probability of failure at any given time may 


depend on how long the component has been functioning satisfactorily. Ex- 
plicitly we have 


component fails during time interval 
10) Pr 4(¢, ¢ + At) if it functions satisfactorilyr = h(é) At. 
up to time ¢ 


Hence if the hazard h(¢) is specified one can easily write down the probability 
density function (equation 9) which represents the failure mechanism. A useful 
consequence of equation (9), is that if 


t 
11) [ r@ ae =, 
0 

then H(¢) follows an exponential distribution with mean-time-between-failures 
equal to @. Thus it is always possible, when dealing with distributions given by 
equation (9), to transform the time-to-failure to a quantity H(#) which will 
follow an exponential distribution. Hence all formulas based upon the exponential 
assumption can be used for the transformed variate H(t). Table 7 summarizes 
the above discussion in terms of a few concrete examples. 

The concept of the instantaneous failure rate is not new and has long been 
used by actuaries. It is often referred to as the “force of mortality’. In con- 
nection with equipment failures this general concept has been pointed out by 
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Carhart [5] and Davis [6]. The relation to other life distributions has been 
discussed by Birnbaum and Saunders [2], Epstein [7], and Zelen [15]. 


6.4 Notes on the analysis and closer approximations to the hypotheses 
H,:bc;; = 1, Hs: a,¢;; = 1, He: a,b;c,; = 1 
The theoretical work from which the techniques presented here were derived 
is given in Zelen [16]. The estimates (equation 5) of the parameters m, a; , b; , and 
c;; are maximum likelihood estimates and the tests of significance in Table 1 
were derived using the likelihood ratio. Tables 2 and 3 were calculated by fitting 
a Pearson curve to the first four moments of the statistics involved and then 


referring to Table 42 in Pearson and Hartley [14] for the normalized percentage 
points. 


TABLE 7 
Examples of Different Life Test Distributions 


Type Conditions An 
of on explicit Probability Transfor- 
Failure hit) hit) density function mation H(t) 





Random: Probability independent h(i) =1/0 1 H® =t 
of failure is inde- of t (t > 0) f® = 9 xP — t/ 

pendent of the time 
that the component 
has been function 
ing satisfactorily. 


(exponential) 


Wear out: Probability increasing h(t) = pte-1/a “ = pie™ an at 

of failure is an in- function (¢ >0,p > 1) 6 

creasing function of of t (Weibull) 

the time that the 

component has been 

functioning satis- 

factorily. 1 1 

Early: Probability decreasing h(t) =e-t/o9 |fJO = ao - 6 {et +1 — ert} H(t) = 1 —e= 
of failure is a de- function (¢ > 0) 
creasing function of t 

of the time that the for small 

component has been t 

functioning satis- 

factorily for small 

values of time. 


(extreme value distribution 
for the range t > 0) 


The test for the hypothesis H, : a;b;c,;; is exactly the test first given by Neyman 
and Pearson and commonly referred to as “Bartlett’s test’’, for testing the 
homogeneity of ab variances each with vy = 2 (r — 1) degrees of freedom arising 
from normal populations. Also it can be shown that the tests for the hypotheses 
H,:b,c;; = land H; :a,c;; = 1 are composed of sums of independent Bartlett 
test statistics. 

Much work has been by Box [3] and Hartley [10] for finding closer approxi- 
mations to the distribution of Bartlett’s test statistic. This can be directly 
carried over to obtain closer approximations on the distribution of the test 
statistics A, and A, . Since A, is exactly Bartlett’s test statistic the application 
of these results is direct. 
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Hartley’s work has resulted in tables [14] which, for our application, give the 
.01 and .05 percentage points of A as a function of the auxiliary quantities k 
and c, . These exact percentage points are those which are labeled (a) in the 
above mentioned tables. On the other hand, Box has shown that a better approxi- 
mation than M = A/C can be obtained by taking f,A/f, (b6* — A) to have an F 
distribution with f, and f, degrees of freedom where f, , f. , and b* are suitably 
defined. Table 8 summarizes the explicit values of the needed auxiliary quantities 
k, c, , and f; , fz , *. 


TABLE 8 


Factors for Closer Approximations to the Distributions of 
Ag, As, and Ag 


Auxiliary quantities Auxiliary quantities! 
Statistic for Hartley’s results for Box’s approximations 


Ay k=b fi = a(b — 1) 
_ vb-1 A b+1 
~ 2b(r — 1) + 6b(r — 1) 
a hi b(a — 1) 
a’ — 1 A a+1 
~ Qa(r — 1) . 6a(r — 1) 
ab f, = ab-—1 


(ab)? — 1 ab +1 
~ Qablr — 1) |- ~ 6ab(r — 1) 


1. The quantities f, and b* are defined by 


i= tt? b* = fe 
. Ay ’ 1 — A, + 2/f2 


2. When the guarantee time G is taken to be zero (G = 0), replace all r — 1 
entries by r. 

3. For the analysis of estimated variances each with v degrees of freedom, 
replace r — 1 by v/2, i.e. r — 1 = v/2. 


7. EXPERIMENTAL SAMPLING TO TEST THE ADEQUACY OF THE 
APPROXIMATIONS TO THE UNKNOWN DISTRIBUTIONS 


The various tests of significance outlined in Table 1 all are of the form M = 
A/C. The quantity A was derived using the likelihood ratio criteria and by 
virtue of this, it is asymptotically distributed as a x’ variate with the appropri- 
ate number of degrees of freedom (when the null hypothesis is true). That is, 
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as the number of observations become large, the probability distribution of A 
approaches that of the x’ distribution. However, for applications, one desires 
the distribution of A for small samples. The problem of finding the exact sampling 
distribution of A is difficult and it is for this reason that small sample approxi- 
mations were resorted to. The small sample approximation is to take M = A/C 
to be distributed as a x” variable, where C is an easily calculated constant. 

The problem now arises as to how well these approximations work. To this 
end, several thousand different life test factorial experiments were generated 
by an electronic computer. The experiments were all taken to have a = 3 and 
b = 5. For each of the ab = 15 treatment combinations, r random variables 
following the exponential distribution with 6 = 1 and G = 0 were generated 
on the computer. (In all cases r was taken equal to n.) The values of r used were 
r = 2,4, and 12. For each value of r, 1000 such 3 X 5 factorial experiments were 
generated. With respect to each experiment the appropriate A and M were 
calculated for the two main effect tests and the test on the interaction. These 
correspond to A, and M, for k = 1, 2, 3 in Table 1. Thus for a particular fixed 
value of r, there will be 1,000 values each of A, and M, for k = 1, 2, 3. 

The probability of a x” variate exceeding A, and M, was calculated for these 
thousand such values, i.e., P{x’ > A,}, P{x’ > M,}. Thus for (say) the test on 
main effects H,(a; = 1), we will have a thousand such probabilities for A, and a 
thousand for M, holding r fixed. If A, follows a x’ distribution exactly, then the 
thousand values of the probabilities are samples from a rectangular distribution 
with range (0, 1). That is when the null hypothesis is true, 5% of the experiments 
would have a level of significance between 0 and 5%; 5% of the experiments 
would be significant between the levels 5 and 10 per cent, etc. A similar remark 
holds for the probabilities associated with M, . 

Figure la summarizes the results of the sampling experiment for the main 
effect test with two degrees of freedom. The histograms on the left side refer to 
the probabilities associated with A, (the likelihood ratio test); the histograms 
on the right hand side refer to the probabilities associated with M, (modified 
likelihood ratio test). These histograms have 20 class intervals. Hence if either 
of the quantities A, or M, follow x’ distributions, we would expect 50 cases 
in each class interval. It is clear from these histograms that the modified likeli- 
hood ratio test (M, = A,/C,) is far superior to the likelihood ratio test (A,) 
for small values of r. When r = 12 it appears as if the asymptotic characteristic 
of the likelihood ratio test is beginning to take hold. Figures 1b and le show 
similar histograms for the main effect test with 4 degrees of freedom and for the 
test on interactions. In all cases it is clear that the approximation M = A/C is 
an improvement over the asymptotic distribution alone. 

In many situations we are interested in the behavior of a test in the tails of 
the distribution. If the distribution of M = A/C is close to that of x’ then we 
would expect that 5% of the generated experiments (50 experiments) would be 
significant at the 5% level of x’. Table 9 summarizes the percent found 
significant both for the modified likelihood ratio and the likelihood ratio test. 
Again the modified likelihood ratio test gives good results and seems fully ade- 
quate for practical applications. 
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TABLE 9 
Percent of Experiments Significant at .05 Level of x? (Exponential Parent) 


Main effects Likelihood ratio Modified likelihood ratio 


2af. 

r= 2 13.8% 
r= 4 7.2% 
r= 12 5.7% 


* 
cae 


16.8% 
9.2% 
6.1% 


22% 
fida 
a : 
NOP bo 


23.8% 
9.2% 
6.5% 


8. RoBUSTNESS 
8.1 Robustness of exponentially derived procedures 


One of the most important characteristics of any statistical technique is 
how well it works when some of the underlying assumptions on which it is 
based are not satisfied. Box and Anderson [4] have introduced the term “robust” 
to describe techniques which are “insensitive to changes in extraneous factors 
not under test’’. In practical applications, one usually does not have enough 
data or a priori knowledge to verify the validity of many of the underlying 
assumptions of a technique. Hence robust procedures are much preferred in 
applications. 

One of the key assumptions of the techniques presented in this paper is that 
the probability distribution of failure times be exponential or more generally 
have the probability density function 


h(t) exp af h(x) az} for t>0 


0 for i< 0. 


In practice, however, one may not know A(?) explicitly. Therefore it is important 
to determine the robustness of the procedures associated with an underlying 
exponential distribution, when in reality the underlying distribution is not 
exponential. 

In order to explore the robustness of the techniques described in this paper, 


sampling experiments were made when the underlying distribution of failures 
was 


12) 3 f(t) = 2te""” (Weibull) 
b) f(é) = exp —{e' — ¢ — 1} (Extreme-value with positive range) 


1) = 
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Figure 2 compares the above probability density functions with that of the 
exponential distribution. 

The sampling experiments consisted of generating 1000 3 X 5 factorial ex- 
periments for r = 2 and 12 having underlying Weibull distributions; for the 
extreme value distribution, 1000 3 X 5 factorial experiments were generated 
with r = 2. The results are summarized in figures 3a and 3b. It is clear from 
these results that the distribution of the significance levels are far from that of a 
rectangular distribution. Hence these techniques are not robust and critically 
depend on the parametric assumption of the failure distribution. 


f(t) = exp -(e*-¢-1) 


£(t) = 2t oe 


2 


3.0 3.5 4.0 4.5 


Figure 2—Comparison of probability density functions for exponential, Weibull, and extreme 
value distributions. 
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8.2 Effect of log transformation on robustness 


It is well-known that if the underlying distribution of failure times is ex- 
ponential, then v6,;/0;; is distributed as a chi-square variate with y = 2 (r — 1) 
degrees of freedom, cf. Epstein and Sobel [9]. Then if we apply a logarithmic 
transformation to the 6;; , we have the model 


13) eo = log 6; = log 0; + log (x’/») 
yi; = log m + log a; + log b; + loge,; + log (x’/r). 


It is also well-known that the distribution of log (x’/v) approaches a normal 
distribution as y — o. Hence an analysis on the y;; using normal theory and 
the usual formulas associated with the analysis of variance seems reasonable. 
Note that the model for y,;; is the usual model for a two-way classification, 
except that the “error term” log (x’/v) has an asymptotic normal distribution. 
The paper by Bartlett and Kendall [1] discusses some consequences of such a 
transformation. 
, In order to compare this procedure with the other procedures, 500 3 X 5 
¢ factorial experiments were made for (i) exponential parent distribution with 
r = 2, 12, (ii) Weibull parent distribution (equation 12a) with r = 2, 12, and 
(iii) extreme value parent distribution (equation 12b) with r = 2. The tests 
on the main effects were carried out using the interaction sum of squares as the 
“error term” in the F-test. The test on interactions was made by using the theo- 
retical variance for log (x’/v) and judging significance by critical values of x’ 
with eight degrees of freedom. 
The results of these sampling experiments are summarized in figures 4a, 4b 
and 4c. From these figures it is obvious that the analysis using the logarithmic 


TABLE 10 . 
Percent of Experiments Significant at .05 Level 


Method* 
of Main effects Main effects Interaction 
Parent Analysis (2 df.) (4 df.) (8 d.f.) 


Exponential M.L. 4.4% 4.6% 7.7% 
log 3.6% 4.6% 8.2% 
M.L. 5.4% 5.8% 5.4% 
log 3.2% 5.2% 4.8% 


Weibull M.L. 2.9% 3.3% 3.4% 
log 5.4% 4.4% 10.8% 
M.L. 0.7% 0.3% 0.2% 
log 5.8% 5.4% 3.8% 


Extreme value M.L. 7.2% 2.7% 3.3% 
log 3.0% 5.2% 6.4% 


* M.L. : Modified likelihood ratio (1000 samples) 
log : Logarithmic transformation (500 samples) 
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transformation results in a startling improvement in robustness. ‘The histograms 
for the Weibull and extreme value distributions are much improved. With 
respect to the tail probabilities, Table 10 summarizes the proportion of experi- 
ments significant at the .05 significance level for all the sampling experiments. 


This table points up very clearly the “robustness” of the logarithmic trans- 
formation. 


Therefore from these sampling studies it is concluded that: 

i) procedures for analyzing life-test data based upon the exponential assump- 
tion depend heavily on this assumption and are not robust; 

ii) preliminary sampling experiments indicate that the logarithmic trans- 
formation may be useful for analyzing life-test data when the underlying failure 
distribution is not explicitly known. 


Acknowledgment: The author would like to express his thanks to Miss M. 
Carroll Dannemiller for programming and supervising the calculations summar- 
ized in sections 7 and 8. 
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Using LaGrange multipliers an example is given wherein the optimum level of 
a fitted second order response function is obtained subject to a constraint provided 
by another fitted second order response function. 


INTRODUCTION 


Quite often one is faced with the problem of having in a given system one 
set of conditions as optimum for a dependent variable y, and another set of 
conditions as optimum for a dependent variable ». A common occurence in 
chemical research is one in which the optimum conditions for yield do not coin- 
cide with those for purity and optimum conditions for one results in an un- 
satisfactory level for the other. A compromise must be made. The best compro- 
mise will ususally be that which gives the maximum yield at some acceptable 
purity level. 

With less than 3 independent variables the surface equations can be solved 
simultaneously by plotting the contour lines for each and the best compromise 
selected by inspection. It is conceivable to extend this to the case of 3 inde- 
pendent variables by the use of solid models or a three-dimensional drawing 
but this would be quite tedious. With 4 or more independent variables this is 
obviously out of the question. 


Tue LAGRANGE MULTIPLIERS 


A solution to this problem has long been known but little recognized by re- 
search workers. It consists of the use of Lagrange multipliers for determining 
the extermal values of a function subject to some constraining relations. 

For example, if 


U = f(a, , te, 23 *** Tq) (1) 
is subject to a constraining relationship 
g(t, ,%2 ,%3 *** Z,) = 0 


and extremum of U is found by solving the equations 


Ue 
na” ee 
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aU dg _ 
ax. +h Or. 


aU dp _ 
a2, + es ee 


o(a, » T (Sey **< Be = 0 


where \ is an unknown multiplier called a Lagrangian multiplier. 

Now if U = f(a, , % , 23) is the equation or mathematical model of our re- 
sponse surface for yield and ¢(x,, 2, 23) = 0 is the explicit function of the purity 
surface for say P% purity, solution of the 4 simultaneous equations (3) in the 
4 unknowns, will give us the conditions for the maximum yield obtainable with 
P% purity product. 


An EXAMPLE 


As a simple illustration let us take a case where 2 independent factors de- 
termine an optimum for both yield and purity but optimum conditions for 
yield are not the same as optimum conditions for purity. We can find the best 
compromise by both graphical and Lagrange methods. 

Consider a response surface for yield adequately represented by 


Y = 55.84 + 7.312, + 26.652. — 3.032? — 6.9623 + 2.692,22 (4) 
Y = predicted yield 

x, = an independent factor such as pressure 

2 = an independent factor such as temperature. 


and a response surface for purity as adequately represented by 


P = 85.72 + 21.852, + 8.59r2 — 9.202? — 5.1822 — 6.262,2, (5) 
P = predicted purity 
2, , 2 = same as for equation (4) 


These equations are plotted as contour lines on figure 1, the upper set for 
equation (4) and the lower set for equation (5). Now if we set a constraint of 
say 90% minimum purity, we see that our maximum possible yield is about 
87% and would occur in the vicinity of x, = 1.0 and 2 = 1.5. Substituting 
these values in (4) gives a predicted yield of 88.43 and substituting in (5) gives 
a predicted purity of 90.21. In solving the problem by the Lagrange method 
we have 

at 
02, 
af 
O22 
aP 


Ox, 


ab 


OX, 


= 7.31 — 6.06z, + 2.692, 


= 26.65 — 13.922, + 2.692, 


= 21.85 — 18.40z, — 6.262, 


= 8.59 — 10,367, — 6.262, 





LAGRANGE MULTIPLIERS WITH RESPONSE SURFACES 


lience, equations (3) become 


7.31 — 6.062, + 2.692, + A(21.85 — 18.40z, — 6.262.) = 0 | 
26.65 — 13.9227. + 2.692, + (8.59 — 10.362, — 6.262,) = 0 


(6) 
85.72 + 21.852, + 8.592. — 9.202} — 5.182; — 6.26.2. = “i 


Equations (6) are non-linear as they will be in the general case. Sometimes 
these equations can be solved by simple iteration. These are several numerical 
methods of solving non-linear simultaneous equations and the non-mathe- 
matician reader is referred to standard texts such as Scarborough [1] or Booth 
[2]. 

Equations (6) in our example were solved by programing the method of 
steepest descent as given by Booth [3] for an IBM 650 Computer. Once pro- 
gramed the optimum yield for any given purity can easily be obtained. Table I 
gives three results for our example which have also been plotted on figure 1. 

The value of such results to a cost analysis of a contemplated process is 


TABLE I 
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quite obvious. The dollar value of premium purity can be balanced against the 
operating cost associated with reduced yields. 
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This paper outlines a technique for determining the probability of failure of a 
plant manufacturing a hazardous product. For the failure of the plant to occur, the 
failure of the safety system must precede the failure of the plant’s critical operating 
system in any “turnaround”’ or inspection period; where the failure of the operating 
system would precipitate a hazardous situation. 

The probability of destructive plant failure is calculated on this conditional basis, 
with probability functions derived accordingly. The derived probability functions— 
which take into consideration the reliability of the individual components compris- 
ing the safety and operating systems, along with selected turnaround cycles—are 
used to illustrate differently designed systems, showing the probable risk of plant 
destruction, where the components are given first in series and then in parallel, 

Under the assumption that the exponential theory of failure applies (for illus- 
trative purposes), probabilities are calculated, showing the effect of changing com- 
ponent reliability and turnaround cycles under given conditions of plant design. 


1. INTRODUCTION 


Chemical plants and petroleum refineries operate a variety of conversion 
processes which are hazardous. Extensive destruction has taken place when the 
operation of these plants inadvertently exceeded normally safe conditions. 
Two recent examples of such occurrences were the destruction of a Fluid Hy- 
droformer unit by gaseous detonation, which “eventually destroyed sixty-three 
tanks and approximately 1,270,000 bbls of crude and various products”, [1] 
and an explosion of a unit handling liquid nitrogen, liquid oxygen, acetylene 
and a newly developed gas with an “estimated damage at about one million 
dollars” [2]. 

The occurrence of such events point to the need for a careful approach to 
safety system design, evaluation, and operation in these industrial areas. D. A. 
Yonge, in his article “Safety in tlh Chemical Industry” [3], for example, lists 
the essential factors in safe plant operation as follows: 


' Paper presented at the Life Testing Reliability Session of the Section in Physical and 


Engineering Sciences of the National Meeting of the American Statistical Association, Chicago, 
Illinois, December 1958. 
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(a) Establishing the requisite safety precautions during the design of process 
equipment. 

(b) Similar considerations with regard to the design of the plant in which the 
process is to be performed with the incorporation of a sufficient number 
and size of safety devices to prevent abnormal conditions or compensate 
for the failure of equipment; 

(c) adequate training of operators, so that each fully understands the possible 
hazards of the process and is fully conversant with safe conditions of 
operation and the steps to be taken in emergency; 

(d) regular and thorough inspection of the plant to ensure that component 
deterioration is detected before failure can occur. 


This paper will provide a method for assessing precautions taken in the opera- 
tion of a plant through an analysis of the design of its safety system together 
with its inspection and maintenance programs. The mathematical model de- 
veloped can assist in the proper selection of plant hazard control systems and 
can provide a means for economic evaluation of incremental risk protection. 


2. Division or A PLANT INTO SAFETY AND OPERATING SYSTEMS 


From the standpoint of plant failure, a manufacturing process can be divided 
into two systems which are amenable to probability analysis: (1) the safety 
system comprised of alarms, operators, automatic trips and valves; and (2) the 


OPERATING COMPONENTS 


ALARM SYSTEM f \ AUTOMATIC-TRIP SYSTEM 


oo nee SYSTEM AUTOMATIC-VALVE SYSTEM Y SYSTEM 


So oS 


FAILURE OF THE PLANT 


WILL OCCUR IF THE SAFETY 
SYSTEM IS IN-OPERATIVE 
WHEN ANY ONE OF THE 
OPERATING COMPONENTS 
HAS FAILED 


Figure 1—Schematic diagram of a typical safety system as related to a set of operating 
components, whose failure precipitates a hazardous situation. 
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operating system comprised of components involved in the production process, 
such as turbines, pipelines, motors and pumps. The objective of the analysis 
is to determine the probability that the safety system will fail and not be able 
to forestall a hazardous situation before a failure of the operating system. 

In the event of a failure in the operating system (here called a primary failure), 
both manual and automatic controls provide for safety, each with their respec- 
tive probability of failure. As illustrated in Figure 1, the manual group consists 
of a system of alarms or visual indicators, calling for the response of the operators; 
the automatic group consists of a system of trips or sensing devices, actuating 
a system of valves or breakers which shuts down the plant and thus prevents 
disaster. For any given primary failure there is a “grace-time”, which is the 
relatively small time-interval between the occurence of the primary failure and 
the subsequent destructive occurence, and enables successful preventive action 
by the operator. 

If either of the two safety groups operate successfully, safety is attained. 
Destruction can only occur if both controls are inoperative, when the primary 
failure takes place. The problem, then, is to evaluate the probability of such 
an event for any proposed safety system as a function of the time between plant 
turnarounds that is, between those periods when the plant is shut down for 
complete inspection and repair. It is assumed that the plant is brought back 
to its initial state at each turnaround. 


3. PROBABILITY OF PLANT FAILURE 


In Figure 2 let the time-axis segment (O, T) represent a turnaround period 
of length 7’. Assuming that the safety system fails at ¢ < T (i.e., during the 
differential time interval ¢ to ¢ + dt), then a hazardous condition will occur if 
the operating system fails in the remaining interval (¢, T). 

Let ¢,(¢) dt denote the probability of safety system failure in the time interval 
(t, ¢ + dt) and let 2,(¢, T’) denote the probability of operating system failure in 
(t, T). 

Since the events are stochastically independent, the combined probability 
of plant failure, during an elementary time interval, is the product of these two 
probabilities denoted by: 


> 
Prob, <X,<t+at {t <a Ss T} = ¢.(t) dtQ.ct, T) ss ¢.(t) at. | c(t) dt 
; (3.1) 


-— | [ 2 ax | dt = g,,,(t, 7) dt 


SAFETY SYSTEM pe OPERATING SYSTEM 
FAILS FAILS 


Oo t t+dt 


Mesa: TURNROUND PERIOD 


(TIME BETWEEN TURNAROUNDS) 


Ficure 2—Safety system failure in (¢, ¢ + dt), and operating system failure in (t, 7’) in a 
turnaround period (0, 7’). 
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where ¢,,.(¢, 7’) dt now represents the joint probability density element for 
plant failure, given ¢,(¢) and ¢,(x) the probability density functions for the 
safety and operating systems, respectively, (using a dummy variable in the 
density function for operating system failure). 

Denote the probability of plant failure in (O, T’) by I'(T), and integrate (3.1) 
over the whole range of ¢. Then the probability of plant failure during a turn- 
around period of length 7’ is represented by the probability function: 


I(T) = [ " g,.dt, T) dt = [ ; 0.0) | aa ax | dt (3.2) 


By permuting the dummy variables in equation (3.2), and changing the order 
of integration, the problem can be formulated so that the destructive failure 
of the plant will occur if the safety system fails in (O, ¢) and the operating sys- 
tem fails in the differential interval (¢, ¢ + dt) as illustrated in Figure 3. The 
probability of plant failure in (O, 7’) is now 


SAFETY SYSTEM OPERATING SYSTEM 
FAILS FAILS 


t t+dt 
TURNROUND PERIOD 
(TIME BETWEEN gem 


Figure 3—Safety system failure in (0, ¢), and operating system failure in (¢, + dt) in a 
turnaround period (0, T). . 


rr) = [ Pee [ ; o.(0] [ “eh ax | dt (3.3) 


where ¢,,-(¢) dé is the joint probability density element. The decision whether 
equation (3.2) or equation (3.3) is to be used depends upon the complexity of 
the safety and operating systems constituting the plant under consideration 
which in turn suggests that a choice between the two formulae will depend upon 
the complexity of the resulting probability density functions for each of the 
systems. When a complicating factor, such as operator failure, enters into the 
probability function for safety system failure, there is an advantage in using 
equation 3.3. 


4. EVALUATION OF THE PROBABILITY OF SAFETY SYSTEM FAILURE 


To derive the distribution function @,(¢) for the failure of the safety system 
in the time interval (0, ¢), the assumption is made that the components com- 
prising the safety system are stochastically independent. It should be pointed 
out that we are interested only in the “first failure’ of these systems, since it 
is reasonable to postulate that no mechanical component can rejuventate itself, 
i.e., that it must either be repaired or replaced, to restore it to working order. 

Denote the probability of failure in (0, ¢) of the manual group (which is 
comprised of alarms and operators) by W4,o(#) and of the automatic group 
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(which is comprised of automatic trips and valves) by W,,y(¢). In accordance 
with the plant design, it will be required that in each system only one component 
has to operate, in order to forestall a failure of that system in the interval (0, #). 
Therefore, 


Viol) = 1— {fl — P.(][1 — Po)]}, (4.1) 
and 


Vr) = 1 — {fl — Pr) — Pr(d]}. (4.2) 


where 


Na 


P,() = I] Pail). (4.3) 


and similarly, one has for the operator, trip and valve systems: 


Py) = I] oily), (4.4) 


Pot) = TT pri(d, (4.5) 


Pt) = TT ped, (4.6) 


given that pu;(t), pri(t), and py,(t) are distribution functions, representing 
the failure of the ith component in each of the systems comprising the manual 
and automatic groups, respectively. 

Since for the failure of the safety system both manual and automatic systems 
must fail, and since these systems are independent, the probability of failure 
of the safety system in (0, #) is given by 


@,(t) = Prob {0 < a, < t} = Wapo(O¥r,r(d) 
= Prob {A or O Fail} Prob {7 or V Fail} (4.7) 
= Prob {AT or AV or OT or OV Fail in (0, )}, 


showing that the failure of the safety system will occur, if at least one of the 
failure events {AT}, {AV}, {OT}, or {OV} occur in (0, #), a process described 
schematically by Figure 4. 

Equation (4.7), equivalent to the following form, provides the distribution 
function for the probability of failure of the entire safety system in (0, /). 


®(t) = {Kq) + [1 — K@)]PsQO} {Pr + Pro) — PrP} = (4.8) 


where K(y) represents a constant denoting the probability of operator failure, 
once P,(y) is determined. 

The probability density function is readily obtained by differentiation of 
equation (4.8): 


d/dt[®,(t)] = ¢.(d). (4.9) 
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a. 

| MANUAL | av Olen 

YSTEM | | SYSTEM 

y 8 
ae — 


EVENTS CAUSING 
SAFETY SYSTEM FAILURE 


Fiaure 4—Events leading to the failure of the safety system. 


The probability of failure of the safety system in (¢, ¢ ++ dé) is now given by 


g(t) dt = “ (Ky) + [1 — KQ))PaO} {Pr + Pr) — PrOPr()}] (4.10) 


5. EVALUATION OF THE PROBABILITY OF FAILURE OF THE OPERATING SYSTEM 


To derive the distribution and probability density functions for the failure 
of the operating system, the assumption is made that the components com- 
prising the operating system are stochastically independent. In addition in 
accordance with the design of the plant, the failure of at least one operating 
component is necessary for the plant to be placed in a hazardous condition. 

Denote the probability of failure in (0, a of the operating system by the 
distribution function 


(i) = Prob {0 <a, < t} =1— I — Dei(t)] (5.1) 


where p,;(t) is the probability of failure in (0, ¢) of the jth operating component. 
The function ®,(¢) is now the distribution function for operating system failure 
in (0, #); hence follows the probability density function 


g(t) = d/di[®.(t)] (5.2) 
and the corresponding probability of failure in (¢, ¢ + dt) of the operating system: 


g(t) dt = at = I aii pa(ol} (5.3) 


6. PROBABILITY OF PLANT FAILURE IN ANY TURNAROUND PERIOD IN 
RELATION TO THE LIFE OF THE PLANT 


Since the probability of failure must be related to the “life of the plant’, 
the probability of plant failure in any turnaround period of length T is evaluated 
as follows: 
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Let n = L/T be the number of turnaround periods during a stipulated plant 


life L. Using the probability T (7) of plant failure for a specified turnaround 
period of length T let 


p= NZ), @=i— i) (6.1) 


respectively, be the probabilities of failure and survival of the plant in (O, T). 


Then the probability of plant failure in any turnaround period of length T can 
be written 


ra, L) =1—- fi - rey, (6.2) 


where 7 is the number of turnarounds in (O, L). For small I'(7), this will be 
approximately equal to nI'(T). 
Since it can also be shown that the probability of the first failure of the plant 
can be expressed as an exponential function, then 
1—T(m, L) =e” (6.3) 


which is the first term of a Poisson series, and in the present problem gives 
the probability of no failure in (O, L), with a failure rate nI'(T). Hence, the 
successive terms for expected r non-failures in (O, L) will be 


Prob {x <r} =e” a , p=ni(T). (6.4) 
1=0 . 


7. Lire oF A PLANT FOR SPECIFIED Risk LEVELS AND 
TURNAROUND CYCLES 


Using equation (6.2), the probability that a plant will survive n turnaround 
periods of length T is given by 


q=1-Ta,D =f{i —ry. (7.1) 


If the probability of survival q and the turnaround period T are specified, 
the life of the plant L = nT can be determined. From equation (7.1), we have 


i pla 
~ dog [1 — T(7)}’ 


so that the life of the plant is given by 


es Wr) < 1. (7.2) 


ad ae : 
L = le ft — TO)’ 0< I(T) <1. (7.3) 


The value L in (7.3) is the life of the plant for which the survival probability 
is q. 


8. PROBABILITY OF PLANT FAILURE, EXPONENTIAL THEORY OF FAILURE 


Although many probability distributions can be used to characterize the 
systems and their respective components in the functions derived above, (4, 
5, 6] the exponential distribution [7] will be used here. In the calculations which 
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follow the integrals were evaluated by Simpson’s 3-point rule, using twenty- 
six intervals. All computations were done on a 205 Datatron Computer. 

Under the assumption that the exponential thedry of failure applies, the 
probability of plant failure may be related to the turnaround, 7’, using the 
density function for the pertinent component: 


ast ee s : 
ee , t>0, a =1/u:, us >0 
0, t<0, 


and the distribution function 


(8.1) 


@,(t) = Prob {a, < t} — ’ t>0, a; = 1/u; , ui > 0, (8.2) 
0, t <0, 

where x; is a continuous random variable denoting the time of failure of the 

ith component, and p; is the mean-time to failure (¢ = A, T, V). 

Using equation (3.3) to determine the probability of plant failure in (O, T), 
we obtain the appropriate joint probability density element ¢, ,.(t)d? by deducing 
the probability differential ¢,(¢) dt for the operating system failure in (¢, ¢ + df) 
and the distribution function #,(¢) for safety system failure in (O, ¢). Thus, 


(i) = {K(y) + [1 — K@)]0 — e***)"} {1 — e87)"* 
+a —e "yh" — (dey — er e")*7} 


is the distribution function for safety system failure in (O, #). Assuming that 
all safety components, comprising a particular system, average the same mean- 
time to failure (i.e., all alarms in the alarm system average the same mean-time 
to failure), equations (4.3), (4.5) and (4.6) can be written 


P(t) = [pO] = l [ “ye * a)" =fi—esy, 


t= A,T,V. 


Using equation (5.1), viz., 


(8.3) 


a) =1- I] - pol 


and substituting 


t 
Dei(t) = [ ane“ dt ie 1 io os 
0 


one finds 
Ne 


@,(i) =1— [Jee =1-e"' 


i 


Hence, follows the probability density function 


el) = d/at [0.0] = (aq) exp { -1> a} 


i 








D) 


6) 
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The probability of plant failure in a turnaround period of length 7 can now 
be determined by substituting the expressions for @,(¢) and ¢.(¢), given in equa- 
tions (8.3) and (8.6) respectively, into equation (3.3). Since cognizance must 
be taken of the “grace-time’’, the turnaround period T should be decreased by 
the increment of time represented by the grace-time y. Hence, 

T=4 
rit, = [ ae*[{Ko) + (d — K@)a — 4} 
, (8.7) 
{(1 ae ae! ti + (i ee eo — (1 bad era ae errr dt 


where 


Ne 
a, = DL a; > N,; is an integer >1, a 20, is A,T, V. 


The period of time .in which a destruction may occur is actually T — y. It 
is readily seen that, if the time period 7 between turnarounds is equal to the 
grace-period , then the period in which destruction may possibly occur is 
equal to T — y = 0, and the integral representing the probability of plant 
failure ['(7’) in (0, T — y) vanishes. Thus, as long asy > T’, the plant can never 
fail. Furthermore since I'(7’) may be evaluated for differently designed safety 
systems involving various numbers of safety components, we attach the sub- 
script s to equation (8.7) to denote the total number of safety components 
comprising the selected safety system. 

Given (8.7) for I',(7, y), we may wish to evaluate the probability of destruc- 
tive plant failure when the safety system comprises only the operator or auto- 
matic groups of the safety system. Thus the safety of the plant is no longer 
protected by a parallel safety system in terms of manual and automatic safety 
groups. Rather than rewrite the basic model, we can obtain expressions for 
these situations by specifying different conditions on the parameters u and a 
and K(7). 

For example, suppose the high cost, of an automatic tripping device, when 
compared with the cost of the plant, precludes the installation of an automatic 
safety system, to supplement or parallel the manual system. Under those cir- 
cumstances, that is, eliminating the automatic system, but retaining the given 
operator-alarm system, we should specify that ur = uy = 0 or, alternatively, 
that ay = ay are very large numbers in the model. Thus 


$() = K(y) + [1 — K@)](l —e**)", since Vr,/y(t) = 1. 


In the same way, the operator-alarm system is eliminated by specifying that 
K(y) = land p, = 0. 

On the other hand, suppose that some specific component within a particular 
safety system should be eliminated from the evaluation. Then we can specify 
B; = ©,a; = Oor, for the operator, K(y) = 0. Thus, any combination of safety 
components or groups can be eliminated from the original model by specifying 
the proper values of the parameters ( or 0) for the pertinent systems. Examples 
of typical modified equations are summarized in Table 1. 
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9. ExaMpLE Usine A SIMPLIFIED PLANT 
As an illustration, consider a simplified manufacturing plant. An evaluation 
of the operating system indicates that two components can be considered as 
“critical” in that their failure will cause extensive plant destruction. To safe- 


guard against this, provision is made for a manual safety system, comprised of 
pertinent alarms and operators. For the sake of simplicity, the operators are 


TRIP-VALVE SYSTEM IS 
CONSIDERED COMPLETELY 
UNRELIABLE AND IS OMITTED 


FROM THE MODEL 


PROBABILITY OF PLANT FAILURE 


Figure 5—Schematic diagram of simplified plant comprised of two critical operating com- 
ponents and a manual safety system. 


considered to be completely reliable, i.e., the probability that an operator will 
fail to take proper action in the stipulated “grace-time”’ is zero. In addition, 


because of high cost an automatic safety system is not installed. Then for the 
purposes of this example the automatic system is considered completely “‘un- 
reliable” and is eliminated from the model. 

In essence, this example illustrates the effect of changing the number of 
alarms in a parallel alarm system, where the governing factors are the average 
failure rate of the components and related turnaround cycles. The problem 
of switching circuits, with the property of shunting an alarm into the initial 
alarm system when any one of them fails, is not being considered here, because 
all components are on stream, and all will operate when an operating system 
failure occurs. 

Following the examples given in Table 1, the probability of plant failure in 
a turnaround period 7, can be determined for the following set of conditions: 


K(y) = wr = py = 0 
a, and ae are finite, 


29 


Vy,(T, 7) = [ ae [1 — e7*4*]"4 dt (9.1) 


Using the numerical data given in Table 2, the probabilities of destructive 
plant failure are graphed in Figure 6, for safety systems containing as many 
as seven parallel alarms, each having two to 150 turnaround cycles during the 
stipulated plant life of fifteen years. 
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TABLE 2 


Numerical Data for the Simplified Plant 
Critical Operating Components 
Ci: Feed Pump 
C.: Feed Pump Motor wee = 10 years 
Safety Components: 
A;: Alarms (i = 1, 2 -::, 7) 
(same mean-time to failure, u,4 , for each of the alarms 
O: Operator 
Stipulated Plant Life: 
L = 15 years 


T(n, 15) 


IN 15 YEARS 


Ma=2 YEARS 
_Me,= 10 YEARS 
Hc,= 10 YEARS 
K(y)=0 


PROBABILITY OF PLANT FAILURE 


20 40 60 80 100 120, 140 
NUMBER OF TURNAROUNDS IN 15 YEARS, n 


Figure 6—Probability of plant failure for simplified plant. Assuming that the underlying 
probability is exponential. 


A first inspection of Figure 6 might lead to the decision to design the safety 
system for a minimum probability of destructive plant failure. Obviously, this 
can only be done at a high cost. Certainly, lower probabilities of destructive 
plant failure can be obtained by inspecting the safety system more often, adding 
more components to the safety system, or make the operating system more 
reliable. However, practical constraints of cost, size, and weight prevent us 
from using an arbitrarily large number of duplicate components [8], and there 
is only a relatively small gain in reliability in a parallel system through the 
addition of more than two components [9, 10, 11]. Thus a loss function must 
be used to determine the least cost combination of turnaround cycles and safety 
devices, while providing a satisfactory degree of surety that the failure of the 
plant will not result in a destructive occurrence. The construction of such a 
loss function will be the subject of a later paper. 


10. ProBABILITY OF DESTRUCTIVE PLANT FAILURE RELATED TO THE 
RELIABILITY OF THE OPERATING SYSTEM FOR THE SIMPLIFIED PLANT 


Frequently, a safety system is so sensitive and reliable that alarms are continu- 
ally being heard throughout the plant everytime that something goes wrong 
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with the operating system. In fact, the safety system is such that it must be 
frequently jammed to prevent its operation. This means that the operating 
system is so unreliable, relative to the safety system that the safety system is 
shutting down the plant with a resulting high loss of production. All operating 
failures, of course, are non-destructive as long as the safety system is not tam- 
pered with. 

Figure 7 is a graph of destructive plant failure relative to the reliability of 
the operating system (reliability measured in terms of the reciprocal mean- 


REGION WHERE PROBABILITY OF PLANT 
FAILURE IS DECREASING, BECAUSE OF 
POOR RELIABILITY OF OPERATING SYSTEM; 
INDICATING THAT OPERATING SYSTEM WILL 
FAIL BEFORE SAFETY SYSTEM 


o 
Co 


IN TURNAROUND OF TWO YEARS 
2 
> 


PROBABILITY OF PLANT FAILURE 


3.0 4.0 5.0 60 de 
0.5 0.33 0.25 0.20 0.166 He, 
RELIABILITY OF OPERATING SYSTEM 


Ficure 7—Graph of probability of plant failure, showing region of non-destructive plant 
failure due to unreliable operating system, resulting in high incidence of plant shutdowns, 
without hazard. Reliability of safety system constant, ag = 2.0. 


time to failure), for a fixed safety system and constant turnaround cycle of 
two years. Reading to the right on the horizontal axis, shows decreasing reli- 
ability of the operating system. The graph shows a maximum probability of 
destructive plant failure at the value a, = 1.2. Between a, = 1.2 and a, = 6.0 
is the region of increasing non-destructive plant failure where the probability 
of operating system failure is increasing up to the point where the plant is 
continually being shut down, and the reliability of the safety system is such 
that it is continually preventing disaster. The further we read along the hori- 
zontal axis, the lower the probability of plant disaster, but the higher the probable 
incidence of involuntary plant shut-downs. 

Thus, to obtain a low probability of destructive plant failure, we must be 
careful to see that such value is obtained with operating system reliabilities 
occurring to the left of the a.-value which marks the maximum probability 
of destructive plant failure. If system failure rates are blindly used for calcu- 
lating the probability of destructive plant failure, a resulting low probability 
for such failure might prove misleading if it occurs in the non-destructive region, 
for this will mean high loss of production. Thus, we should strive for low proba- 
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T,(2;X) 
1.0 


ENVELOPE WHICH CONTAINS 
REGION OF NON-DESTRUCTIVE 


PLANT FAILURE 


o o o 
oS o oo 


TURNAROUND OF TWO YEARS 
i 
Nu 


PROBABILITY OF PLANT FAILURE IN 


12. 13 
RELIABILITY OF OPERATING SYSTEM, X=At2a,, A=2aq 


Ficure 8—Graph of probability of plant failure, showing envelope containing region of non- 
destructive plant failure due to changing reliability of the operating and safety systems. 


bilities of plant failure which occur in the shaded region illustrated by the graph 
of Figure 7. 

Figure 8 is a graph of destructive plant failure (given a constant turnaround 
cycle of two years) relative to varying reliabilities of both safety and operating 
systems, showing the envelope which contain the region of non-destructive 
plant failure due to increasingly unreliable operating systems. 

If we now augment the alarm system with an automatic system, comprised 
of a single trip and valve, the resulting envelopes containing the region of 
non-destructive failure plant shifts accordingly. These envelopes are shown in 
Figure 9. 


r [Na. Nr. Ny] (2#*) 
10 


=2 


NOTE: ay=ay=0, IN A, FOR SYSTEMS 
WITH Ny=Ny=0 


IN A TURNAROUND OF TWO YEARS, T 


2.0 3.0 4.0 5.0 6.0 7.0 8.0 
X=At+2ac, A=2(aatattay) ; 


PROBABILITY OF HAZARDOUS PLANT FAILURE 


Figure 9—Graph of probability of hazardous plant failure in (0,2), showing envelopes for 
differently designed safety systems. 
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11. ConcLUsION 


This paper has outlined a technique for determining a finite probability of 
destructive failure of a plant producing a hazardous product. Such failure is 
a function of the reliability of two independent and parallel safety systems, 
and specified critical operating components, whose failure which would cause 
extensive plant destruction. 
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